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This  report  presents  and  compares  methods  for  computing  the 
gravity  vector  outside  the  earth.  The  gravity  vector  is  conveniently 
split  up  into  a  normal  part,  the  vector  of  normal  gravity,  and  an 
anomalous  part,  the  vector  of  gravity  disturbance. 

Part  I  gives  the  theoretical  foundations  and  a  practical  method^ 
for  the  computation  of  the  vectors  of  normal  gravity  and  gravitation. 

The  method  is  adapted  for  electronic  computation.  It  is  illustrated 
by  a  numerical  example 

Then  the  formulas  for  the  computation  of  the  gravity  disturbance 
vector  from  free-air  gravity  anomalies  Ag  at  the  surface  of  the  earth 
are  developed,  (a)  for  the  direct  method,  which  uses  Ag  only,  and 
(b)  for  the  coating  method,  which  requires  Ag  and  the  geoid  heights  N 
but  involves  simpler  formulas. 

The  components  of  the  gravity  vector  are  first  computed  in  a  local 
coordinate  system  and  then  transformed  to  the  geocentric  world  system 
by  a  spatial  rotation.  If  the  gravity  vector  is  required  along  a  rocket 
trajectory,  then  also  the  components  along  and  normal  to  the  trajectory 
can  be  computed  by  a  spatial  rotation, 

tV 

Part  II  is  concerned  with  accuracy  studies_,in  order  to  determine 
the  best  practical  procedure  for  the  computation  of  the  gravity  disturbances. 

The  standard  errors  of  the  components  of  the  gravity  disturbances, 
due  to  the  interpolation  of  the  gravity  material,  are,  approximately, 
inversely  proportional  to  the  elevation  and  very  small  for  high  altitudes, 
provided  there  is  uniform  coverage  by  gravity  stations. 

The  influence  of  the  distant  zones  is  considered  in  some  detail. 

This  influence  decreases  very  slowly  beyond  a  certain  radius,  so  that  it 
is  impracifical  to  go  farther  than  about  30°  in  "the  direct  method  and  20° 
in  the  coating  method  (with  an  error  of  about  ±5  mgal  in  both  cases), 
unless  the  integration  is  extended  over  the  whole  earth,  which  is  necessary 
for  higher  accuracy. 

There  is  another  method,  for  which  the  influence  of  the  distant 
zones  is  completely  negligible  and  which  furthermore  is  the  simptest. 
the  upward  continuation  of  the  surface  disturbances.  If  presupposes. 
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however,  the  deflections  of  the  vertical  5  and  T]  for  the  horizontal 
components,  and  Ag  and  N  for  the  vertical  component. 

For  practical  use  the  following  methods  are  proposed:  if  only 
Ag  is  given,  the  direct  method;  if  Ag  and  N  are  given,  the 
coating  method  for  the  horizontal  components  and  the  upward  continuation 
for  the  vertical  component;  if  Ag  ,  5  ,  T|  are  given,  the  upward 

continuation  for  all  three  components. 

Finally,  a  detailed  practical  computation  procedure  is  described 
for  the  practically  most  important  case  that  Ag  and  N  are  given. 


II 


TABLE  OF  CONTENTS 


Part  I 

by  R.  A,  Hirvonen 

Introduction  .  1 

1.  Conversion  of  Coordinates  .  1 

2.  The  Potential  of  Normal  Gravitation  .  9 

3.  The  Components  of  Normal  Gravity . 10 

4.  Summary  of  Formulas  for  Normal  Gravity . 13 

5.  Gravity  Disturbances  .  l8 

Part  II 
by  H.  Moritz 

6.  Influence  of  an  Inaccurate  Gravity  Material  ...  27 

7.  Influence  of  the  Remote  Zones . 38 

8.  Practical  Computation  of  the  Gravity  Disturbances  .64 

References . 87 


III 


PART  I 


by 

R.  A.  Hirvonen 

INTRODUCTION 

The  purpose  of  this  report  is  to  develop  a  method  of  computation 
of  the  gravity  vector  for  a  great  number  of  points  at  high  altitudes, 
e.g.  for  trajectories  of  moving  bodies,  when  the  geocentric  coordinates 
are  given.  Formulas  wdll  be  summarized  in  a  form  suitable  for  auto¬ 
matic  high  speed  computers.  Nmerical  examples  illustrate  the  method. 

First,  the  components  of  the  normal  gravitation  (attraction  of  a 
reference  ellipsoid)  will  be  computed.  Those  of  the  normal  gravity 
can  then  be  obtained  by  subtraction  of  the  components  of  the  centri¬ 
fugal  force.  Finally,  the  deviations  of  the  actual  components  from 
these  normal  values  will  be  evaluated  on  the  basis  of  the  free  air 
anomalies  of  the  gravity  observed  on  the  physical  surface  of  the  earth. 
The  accuracy  of  the  method  is  studied  in  Part  II. 

The  method  is  based  on  the  formulas  first  published  in  [Hirvonen, 
1959],  hut  several  practical  improvements  have  been  introduced. 

1.  CONVERSION  OF  THE  COORDINATES 


It  is  supposed  that  the  Cartesian  geocentric  coordinates  x,y,z 


are  given  for  each  point  P  at  which  the  gravity  will  be  computed . 

The  z-axis  is  the  axis  of  the  earth's  rotation^  the  x-axis  has  the 
longitude  0°  (Greenwich)  and  the  y-axis  has  the  longitude  90°E.  These 
coordinates  must  be  converted  into  three  other  systems. 

I.  Geocentric  coordinates,  defined  by 
X  =  r  cos  Y  cos  X 
(l)  y  =  r  cos  f  sin  X 
z  =  r  sin  f 
and  called 

r  =  radius  vector 
f  =  geocentric  latitude 
X  =  geographic  longitude 

When  x,y,z  are  given,  r,Y, X  can  easily  be  computed: 


(1.1) 

P"  = 

=  x= 

+  j 

(1.2) 

r®  = 

=  P" 

+  z 

(1.3) 

cos 

X  = 

X 

p 

(I.U) 

sin 

X  = 

y 

p 

(1.5) 

cos 

y  = 

p 

r 

(1.6) 

sin 

y  = 

z 

r 

For  the  computation  of  the  normal  gravitation,  we  shall  use  the 


function 


2 


2  2 

(1.7)  F  =  cos  2  T  =  ^  ~g^ 

II.  Geographic  coordinates,  defined  by 
p  =  (N  +  H)  cos  cp 

(II) 

z  =  (N  +  H  -  e^N)  sin  cp 

(1.8)  N  =  I 

w 

(1.9)  1  -  e^  sin^  cp 
a  a®-b^ 

and  called 


cp  =  geographic  latitude 
H  =  geographic  height 

The  geographic  coordinates  refer  to  a  reference  ellipsoid  vith  equa¬ 
torial  radius  a  and  polar  radius  b  ,  and  they  are  used  for  the 
computation  of  triangulations,  in  spite  of  the  fact  that  the  "orienta¬ 
tion"  of  the  ellipsoid  (or  the  "data"  of  the  triangulation  system)  can¬ 
not  be  determined  without  world-wide  geodetic  operations.  The  names 
suggested  here  are,  however,  not  yet  in  a  general  use. 

The  normal  gravitation  is  the  attraction  of  the  reference  ellip¬ 
soid  which  is  supposed  to  have  the  same  mass  as  the  actual  earth.  The 
masses  inside  the  rotating  ellipsoid  are  supposed  to  be  distributed 
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in  such  a  way  that  the  combined  potential  of  the  attraction  and  of  ' 
the  centrifugal  force  is  constant  at  the  outer  surface  and  very  close 
to  the  actual  potential  of  the  gravity  at  the  mean  sea  level. 

At  higher  altitudes,  the  geographic  coordinates  are  less  practical 
for  the  exact  computation  of  the  geodetic  quantities.  Especially, 
we  cannot  use  them  for  the  computation  of  the  normal  gravity.  However, 
we  should  use  them  for  the  evaluation  of  the  deviations  of  the  actual 
gravity  from  the  normal  gravity  on  the  basis  of  the  anomalies  observed 
on  the  ground. 

The  computation  of  cp  and  H  from  x,y,z  is  rather  complicated. 


The  difference  of  these  equations  gives 
z  W 

(1.12)  tan  cc  =  —  +  e^  —  sin  cp 

P  P 

This  equation  will  be  solved  with  respect  to  cp  by  successive  approxi¬ 
mations.  The  first  approximation  is  usually  obtained  by  H  =  0  in  which 
case  we  have 

(1.13)  tan  CTi  =  ^ 

l~e  p 
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For  the  higher  altitudes,  it  is  slightly  better  to  start  from 


ular  purpose.  However,  we  have  included  the  computation  of  the  correct 
value  into  our  program  because  it  may  be  useful  for  other  purposes. 

Note  that  cp  is  not  the  direction  of  the  normal  gravity  at  the 
elevated  point  P  but  the  direction  of  the  normal  of  the  reference 
ellipsoid.  The  former  direction  can  be  computed  only  when  the  compo¬ 
nents  of  the  normal  gravity  are  known: 

(1.20)  tan  cp"  = 

y 

p 

We  shall  call  cp"  the  geodetic  latitude,  bearing  in  mind  the  definition 
§  -  component  of  the  deflection  of  the  vertical  is  the  difference  of 
the  astronomical  latitude  cp'  and  geodetic  latitude  cp"  at  the 
elevated  point  P  . 

The  difference  cp"-  cp  is  caused  by  the  curvature  of  the  normal 
plumb  line.  In  [Hirvonen,  i960],  a  series  is  given  for  the  computation 
of  cp"-  cp 

cp"  -  cp  =  sin  2  cp  (0",  170293  H 

(1.21)  +  0  ,001103  H  cos  2  cp 

+  0  .000034  H^) 

(H  in  kilometers).  In  the  present  problem,  however,  we  cannot  use 
this  formula  for  the  computation  of  cp  because  H  still  is  unknown. 
Therefore,  we  have  to  compute  cp  by  successive  approximations  as 
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described  above.  When  the  final  approximation  has  been  founds  H  can 
be  obtained  from  (l.lO). 

For  latitudes  higher  than  45°,  our  program  should  be  based  on  an 
alternative  set  of  formulas.  Instead  of  (l.l4)  through  (I.I9),  compute 


(1.22)  cot  cpi  - 


(1.24)  cot  (fg  = 


p  e  a  cos  cpi 

z  Wi  z 


=  cot  Y  -  5 


(1.25)  cot  cp3  -  cot  cfjg  =  6  •  (gQ-t  qjg.  cot  cpi ) 

COuCPj, 

With  the  final  approximation  of  cp,  compute  H  from  (l.ll), 
III.  Elliptic  coordinates,  defined  by 


p  = -  cos  p 

sin  e 


(III)  z  =  ■  sin  p 


tan  e 


c  =  e  a 

and  called 

e  =  angular  eccentricity 
P  =  reduced  latitude 

At  the  surface  of  the  reference  ellipsoid,  e  has  a  constant  value 
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e  with 
o 

(1.26)  sin  e  =  e 

o 

When  p  and  z  are  given,  we  could  compute  e  and  p  by 
rigorous  formulas 

(1.27)  k®  =  r^  +  c® 

(1.28)  h^  =  k^  -  i^p^ 

-1,2 

(1.29)  sin^  e  =  '■  ~a 

2p‘’ 

(1.30)  cos®  P  = 

For  practical  computations,  however,  we  shall  use  power  series. 
By  the  aid  of 

(1.31)  H®  =  ^ 


we  can  replace  the  formulas  above  by 

(1.32)  k®  =  r®(l  +  K®) 


(1.33) 

11 

r^ 

H 

1 

ro 

®  cos 

2Y  + 

(1.34) 

sin® 

e 

=  K®[1 

-  H® 

sin®  Y 

- 

sin® 

Y  cos  2  Y 

-  h6 

.  2 

sin 

Y  (7-  cos 
4 

2 

1  1 

t— 1  1 

(1.35) 

cos® 

P 

2 

=  COS 

Y  [1 

-  H®  sin® 

Y 

sin®  Y 

cos 

2  Y] 
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Using  the  abbreviating  symbol 


(l.T)  F  =  cos  2  Y 

we  obtain,  after  lengthy  but  easy  computations, 


e 

(1.36) 


H  +  ~  (3F  -  1) 

+  k5(35F®  -  lOF  -  13) 

+  h^(33F^  -  9F®  -  21F  +  —  ) 

128  ^  7 

2.  THE  POTENTIAL  OF  NORMAL  GRAVITATION 


The  Newtonian  attraction  of  the  reference  ellipsoid  is  called  the 
normal  gravitation  of  the  earth.  The  potential  of  this  attraction  can 


be  expressed  in  a  closed  form: 

/  \  uu^a®  ,2  2  \ 

(2,1)  V  =  —  e  +  q(-  -  cos^  P) 


where 


f  :  gravitational  constant, 

M  :  mass  of  the  earth, 

u)  :  rotation  speed  of  the  earth. 

The  aiixiliar  variable 

(2.2)  q  =  ■^  [e  -  3  cot  e  (1  -  €  cot  e)] 

is  a  function  of  e  only.  At  the  surface  of  the  reference  ellipsoid 

q  has  a  constant  value  q  . 

o 
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The  closed  formula  (2.2)  of  q  is  not  suitable  for  practical 


computations,  but  must  be  replaced  by  power  series: 

234527 
(2.3)  "  15  ^  ^  «  -  .  .  . 


or 


/o  k'*  2.3  3.5  5  sin'  c  +  .  . 

(2.4)  q  =  —  sin  c  +  —  sin  c 


Therefore,  it  is  only  logical  that  we  use  power  series  throughout  for 
our  computations.  ?y  the  aid  of  (1.34),  (2.4)  gives 

(a.5)  1  -  ^  ^  -  7)  -  ^  (63J*  -  lOF 

Now  we  have  to  insert  (I.35),  (1.3*5)  and  (2.5)  into  (2.1).  Using 
another  abbreviating  symbol  for  the  constant 


(2.6)  C  =  fr 


2  <«®a*c 


15  fM^ 


we  can  write  the  result  in  form 


V  =  ^  h{i  +1  kMi  -  C)(F  -  i) 

(2,7)  -  10C)(7F®-  2F  -  ^) 

+  (3  -  5C)(11F^-  3F®-  7F  +  ^)} 


3,  THE  COMPONENTS  OF  NORMAL  GRAVITY 


Differentiation  of  (2„7)  and  (I.3I)  with  respect  to  r  gives 
the  component  of  the  normel  gravitation  toward  the  center  of  the  earth: 
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(1  -  C)(3F  -  1) 


^  f M  2  .  ^  ; 

H  [1  +  —  H 


r  =  -  .  =  a 

r  or  c 


1  2 


(3.1) 


+  ^  (  7  -  10C)(35F2  -  lOF 

+  -^  (3  -  5C)(231f''  -  63f2 


13) 

l^TF  +19)} 


Differentiation  of  (2.7)  and  (I.7)  with  respect  to  7  and 


division  by  r  give  the  component  perpendicular  to  p  ■ 

1  5F  FM  4  .  ^  Mf  n 

=  ■  7  ii  =  ^  ^  sxn  Y  cos  Y  {1  -  C 


(3.2) 


+  ^  H2  (7  -  10C)(7F  -  1) 


+  ^  (3  -  5C)(33F^  -  6F  -  7)} 


In  order  to  facilitate  the  numerical  computations,  we  replace 
the  increasing  powers  of  by  those  of 

(3.3)  G  =  ;j 

using  convenient  units.  In  other  words,  we  insert: 


(3.1.)  = 


For  the  international  ellipsoid  with  the  international  normal  formula 

of  gravity  we  obtain,  when  r  is  expressed  in  kilometers  and  F  in 

-2 

gal  =  cm  sec  ; 

(3.5)  Q  ^  3986  3290.49 

r= 
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r  =  1000  G 
r 

-  0.835  888 

+  2.507  664  G^  cos  2Y 

-  0,005  118  G^ 

-  0.003  937  0^  cos  27 
+  0.013  778  G^  cos®  27 
+  0.000  007  G^ 

-  0,000  054  G^  cos  27 

-  0,000  023  G^  cos®  27 

.  a 

+  0,000  085  G^  cos  27 


r 


7 


(3.7) 


sin  7  COE  7  (3.3^3  551  G® 

-  0,0C3  i49  G® 
+  0,022  045  G® 

-  0,000  031  G^ 

-  0,000  026  G^ 

+  0.000  145  G*^ 


cos  27 

cos  27 
cos®  27} 


As  G  is  always  slightly  smaller  than  1,  it  is  easy  to  see 
that  the  accuracy  of  one  milligal  can  he  obtained  without  the  terms 
with  G^ . 

The  components  of  gravitation  in  directions  of  the  Cartesian 
coordinate  axes  are  here  chosen  to  he  positive  towards  the  origin. 
They  can  he  computed  hy  formulas 

(3.8)  r  r  cos  7  -  r  sin  i 

p  r  7 


(3.9)  r  =  r  cos  X 

X  p 
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(3.10)  r  =  r  sin  X 

y  p 

(3.11)  r  =  ^  sin  Y  +  r  cos  Y 

z  r  I 

The  components  of  the  normal  gravity  are  obtained  by  subtraction 
of  those  of  the  centrifugal  force: 

(3.12)  y  =  r  -  uj2p 

P  P 

(3.13)  y  =  r  -  U)®x 

X  X 

(3.14)  y  =  r  -  u)^y 

y  y 

(3.15)  y  =  r 

z  z 

The  total  gravitation  is 

(3.16)  r  =  /r  r  * 

r  I 

and  the  total  gravity 

(3.17)  y  =  /y  ^-hy  ^ 

P  z 

The  direction  of  the  former  is 

p 

(3.18)  §  =  arc  tan 

P 

that  of  the  latter  was  given  in  (l.20), 

4.  SUMMARY  OF  FORMULAS  FOR  NORMAL  GRAVITY 
The  formulas  given  above  will  be  summarized  here  in  form  of  a 
program  for  automatic  computers.  First,  the  constants  are  listed  which 
pertain  to  the  international  ellipsoid  and  to  the  international  normal 
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formula  of  gravity.  Then  the  computations  are  described  in  detail^ 


using  a  numerical  example  as  illustration. 


Input 

Constants 

(1) 

1,0000  0000 

(2) 

e^ 

=  0,0067  2267 

(3) 

1  +  e® 

=  1,0067  2267 

ih) 

a 

=  6378,3880 

(5) 

e«a 

=  14-2,8798 

(6) 

f  M 

=  3986  3290 

(7) 

1000,0000 

(8) 

0,83589 

(9) 

2,50766 

(10) 

0,00512 

(11) 

0,0039^ 

(12) 

0,01378 

(13) 

3,34355 

ilk) 

0,00315 

(15) 

0,02204 

(16) 

0)= 

=  0,000  531  749 

Input 

Variables 

Example 

(IT) 

Identification  number  of  P 

(18) 

X 

521,051 r 

(19) 

y 

4957.4762 

(20) 

z 

4169.1994 

Ik 


Program  of  Computations 


(21) 

(18)  •  (18)  +  (19)  •  (19) 

2484  8065 

(22) 

(20)  .  (20) 

1738  2224 

(23) 

r® 

(21)  +  (22) 

4223 ,0289 

(24) 

cos  2Y 

[(21)  -  (22)]  V  (23) 

0.1767  8877 

(25) 

P 

/(21) 

4984,7833 

(26) 

r 

/(23) 

6498,4836 

(27) 

cos  7 

(25)  -  (26) 

0.7670  6869 

(28) 

sin  7 

(20)  ^  (26) 

0.6415  6496 

(29) 

tan  7 

(20)  (25) 

0.8363  8528 

(30) 

cos  X 

(18)  -  (25) 

0.1045  2846 

(31) 

sin  X 

(19)  -  (25) 

0,,9945  2191 

(32) 

cot  X 

(18)  V  (19) 

0.1051  0423 

(33) 

cosY  cosX  (27)  •  (30) 

0.0802 

(34) 

sinY  cosX  (28)  •  (30) 

0,0670 

(35) 

cosY  sinX  (27)  •  (31) 

0.7629 

(36) 

sinT  sinX  (28)  •  (31) 

0.6381 

(37) 

X 

arc  cot  (32) 

84° 00;  0000 

(38) 

tan 

(3)  •  (29) 

0.8420  0802 

(39) 

tan®  cpi 

(38)  .  (38) 

0.7089  7751 

(40) 

sin®  cpi 

(39)  V  [(1)  +  (39)] 

0.4i48  5479 

(41) 

W® 

(1)  -  (2)  •  (40) 

0.9972  1107 

(42) 

sin  cpi 

/(40) 

0.6440  9222 

(43) 

W 

/(4l) 

0.9986  0456 

(44) 

(5)  •  (42)  V  (25) 

0,0055  4057 

(45) 

S 

(44)  V  (43) 

0,0055  4831 
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(46) 

tan  (ffe 

(29)  +  (45) 

(47) 

(38)  -  (46) 

/  (48) 

(38)  •  [(1)  + 

(49) 

(47)  •  (45)  - 

(50) 

tan  CD 

(46)  -  (49) 

(51) 

tan®  cp 

(50)  •  (50) 

(52) 

sec  9 

y(l)  +  (51) 

(53) 

cos  cp 

(1)  V  (52) 

(54) 

sin  cp 

(50)  •  (53) 

(55) 

9 

arc  tan  (50) 

(56) 

(25)  •  (52) 

(57) 

(4)  -  (43) 

(58) 

H 

(56)  -  (57) 

(59) 

G 

(6)  -  (23) 

(60) 

G2 

(59)  •  (59) 

(61) 

G®F 

(60)  .  (24) 

(62) 

3 

G 

(60)  •  (59) 

(63) 

G^F 

(62)  •  (24) 

(64) 

G^F2 

(63)  •  (24) 

(65) 

(7)  •  (59) 

(66) 

(8)  •  (60) 

(67) 

(9)  •  (61) 

(68) 

(10)  •  (62) 

(69) 

(11)  •  (63) 

(TO) 

(12)  •  (64) 

0.8419 

3359 

0,0000 

7443 

(39)] 

1.4390 

(48) 

0,0000 

0029 

0.8419 

3330 

0 

-Tj 

0 

CO 

CO 

5168 

1.3072 

3054 

0.7649 

7601 

0.6440 

5878 

4o°05 ' 

.7085 

6516.2610 

6387.3011 

128.96 


0.9^39  5022 
0.8910  k2 
0.1575  26 
0,84ii 
0,1487 
0.0263 
943.950  22 
0.744  81 
0.395  02 
0.004  31 
0.000  59 
0.000  36 
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Input 

Variables 

Example 

(71) 

(65)-(66)+(67)-(68)-(69)+(70) 

943.595  90 

(72) 

(13)  •  (60) 

2,979  24 

(73) 

(14)  •  (62) 

0,002  65 

(74) 

(15)  •  (63) 

0.003  28 

(75) 

(72)  -  (73)  +  (74) 

2,979  87 

(76) 

(75)  •  (27)  •  (28) 

1,466  47 

(77) 

^p 

(71)  •  (27)  -  (76)  •  (28) 

722,862  o4 

(78) 

(71)  •  (28)  +  (76)  •  (27) 

606,502  95 

(79) 

(77)  •  (30) 

75.559  66 

(80) 

(77)  -  (31)  - 

718,902  l4 

(81) 

^p 

(77)  -  (16)  •  (25) 

720,211  39 

(82) 

V 

X 

(79)  -  (16)  •  (18) 

75.282  59 

(83) 

V 

y 

(80)  -  (16)  .  (19) 

716,266  01 

(84) 

tan  $ 

(78)  t  (77) 

0,8390  3002 

(85) 

tan  cp" 

(78)  t  (81) 

0.8421  1796 

(86) 

fS 

(71)  •  (71)  +  (76)  •  (76) 

890375.37 

(87) 

a 

r 

(78)  .  (78)  +  (81)  .  (81) 

886550,27 

(88) 

r 

/(86) 

943.59704 

(89) 

y 

/(87) 

941.56799 

(90) 

arc  tan  (84) 

39°59’ .8596 

(91) 

cp" 

arc  tan  (85) 

4o°o6' .0800 

Output 

(17), 

(55),  (37), 

(58),  (53),  (79),  (80),  (78) 

in  full 

(33), 

(34),  (31), 

(35),  (^),  (30),  (28),  (27) 

four  digit  values 
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5.  GRAVITY  DISTURBANCES 

If  the  free  air  anomalies  Ag  are  knovn  everywhere  on  the  physical 
suhface  S  of  the  earth,  the  potential  disturbance  T  at  any  point 
on  or  above  S  can  be  computed.  The  computation  is  very  complicated, 
especially  if  the  point  is  very  close  to  high  and  steep  mountains.  In 
most  cases,  however,  the  generalized  formula  of  Stokes  gives  a  sufficient 
approximation . 

We  shall  use  following  notations: 

P  :  fixed  point  at  which  T  is  wanted; 

M  ;  "moving"  point  at  S  ; 

r  :  distance  from  the  centre  of  the  earth; 


(5.2)  D^  =  1  -  2t  cos  Y  +  t^; 

(5.3)  d.a  =  sin  y  'iCt:  =  cos  cp  dcp  dx 

is  the  element  of  solid  angle,  situated  at  point  M  . 

The  generalized  formula  of  Stokes  reads 

(5-^)  J  ^  [-  +  1  -  3D  -  t  cos  Y  (5+3^n - ^ - )}da 
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The  integration  must  he  carried  out  over  the  entire  surface  S  .  To  a 
very  good  approximation^ 

(5.5) 

where  H  is  the  height  of  P  above  S  .  Usually,  r  can  even  be 

M 

replaced  by  a  constant  R  ,  the  mean  radius  of  the  earth. 

Especially,  if  H  =  0  or  t  =  1,  we  have 
(5  .6)  D  =  2  sin 
and 


(5-7)  To-r;  jAs^,  { 


w  r 

,  ,  cosec  -^  +  1  -  D  sin  — 

2  2 


-  cos  '1'  [5  +  3  (sin  +  sin®  ■^)]}da 

2  2 

The  quantity 

(5.8)  c  =  ; 

is  called  the  height  anomaly  and 

is  approximately  the  elevation  of  the  geoid  above  the  reference  ellipsoid. 

Because  C  is  one  of  the  most  important  quantities  of  geodesy,  we 
S 

may  often  assume  that  the  values  of  it  have  already  been  computed. 


Then  we  can  compute  the  "coating"  function  of  Helmert: 

3  Cq 

(5.10)  n  =  Ag  +  I  -f 
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This  function  has  the  advantage  that  the  long  formula  (5.^)  for  the 


computation  of  T 
(5.11) 

P  STTTp 


at  higher  altitudes  can  he  replaced  hy  a  shorter  one 


In  our  present  problem^  we  have  to  compute  the  components  of  the 
gravity  disturbance: 


(5.12) 

(5.13) 
(5.li^) 

If 

written 


n  3r 


cos  q:  5T 
r  a'f 


sin  a  St 


we  take  T  from  the  long  formula  (5.^),  the  result  can  be 
in  form 


(5.15)  \  =  ^  da 

(5.16)  6^  =  ^  J Ag  Fs  cos  a  da 

(5.17)  6^  =  I;;  J Ag  Fa  sin  a  da 


where 

(5.18)  Fi  =  i  t^  ^  +  1  -  6d 

-  tcosY  (i3.6£n^ii^|^^)} 

(5.19)  Fa  =  t®  sin  Y  {p  +  |  -  4 


^  /F-d+'t  cos  Y 
2  D  sin®Y 


Zn 


IM-l-t  COS  Y 


)} 
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If  we  use  the  short  formula  (5.11),  we  obtain 

(5.20)  5^  =  I;  J  1*  fiaa 

(5.21)  6  =  ■^  r  i’s  cos  a  do 

^  m  2tt  'J  ^ 

(5.22)  5 ,  =  —  r  M-  fo  sin  OL  da 

'  2tt  ^ 

where 

(5.23)  f^  =  ^  (1  -  t  cos  Y) 

t"^ 

(5.24)  i’s  =  p-  sin  Y 

Table  5.1  shows  some  values  of  and  fi,  Table  5.2  those  of  Fg 
and  fa .  We  see  that  for  small  values  of  Y  ,  the  functions  F  and  f 
are  almost  equal. 

In  practical  computations,  we  must  first  find  the  values  of  Y 
and  a  .  If  the  earth  is  considered  as  a  sphere,  we  have 

(5.25)  cos  Y  =  sin  m  sin  cp  +  cos  cp  cos  cp,,  cos  (x  -  \) 

P  M  P  M  M  P 

(5.26)  sin  Y  cos  a  =  cos  cp^  sin  cp^^  -  sin  cp^  cos  cp^  cos  (4^-  ^p) 

(5.27) -  sin  Y  sin  a  =  cos  cp^^  sin  (4j^-  ^p) 

For  small  values  of  Y  ,  we  can  use  the  approximate  formulas: 

(5.28)  m  =  cpj^  -  cpp 
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Tatle  5.1 


H  =  2  km 

H  =  8 

Y 

0' 

10l6xxxx 

636XXX 

1 

4009xxx 

588xxx 

3 

395  3xx 

3522XX 

10 

132IX 

4002x 

30 

710 

2062 

1° 

182 

356 

3  1 

45 

51 

10 

12.1 

12,2 

30 

-  3.7 

-  3.7 

60 

-  3.3 

-  3.3 

90 

-  2.3 

-  2.3 

120 

0.4 

0.4 

150 

3.0 

3.0 

l80 

4.1 

4.1 

O' 

10l6xxxx 

635  XXX 

1 

4005 XXX 

587XXX 

3 

3937XX 

3512XX 

10 

1267X 

395 6x 

30 

528 

1883 

1° 

88 

262 

3 

12 

18 

10 

2.9 

3.1 

30 

1.0 

1.0 

60 

0.5 

0.5 

90 

0.4 

0.4 

120 

0.3 

0.3 

150 

0.3 

0.3 

180 

0.2 

0.2 

H  -  32 

CO 

OJ 

1 — 1 

il 

w 

H  =  512 

400xx 

2559 

174 

398xx 

2558 

174 

3824x 

2552 

174 

259^x 

2481 

173 

5085 

1980 

171 

947 

1120 

162 

76 

i48 

105 

12.8 

14.9 

18.8 

-  3.6 

-  3.^ 

-  2,7 

-  3.3 

-  3.1 

-  2.4 

-  2.3 

-  2.2 

-  1.8 

0.4 

0.4 

0.3 

3.0 

CO 

CM 

2.3 

4.0 

3.8 

3.1 

397XX 

2480 

155 

395  XX 

2480 

155 

3795  X 

2472 

155 

2568X 

2403 

154 

4928 

1908 

152 

858 

1061 

l44 

44 

118 

90 

3.8 

6.3 

12.2 

1.0 

1.1 

1.3 

0.5 

0.5 

0.5 

0.4 

0.3 

0.3 

0.3 

0.3 

0.3 

0.3 

0.3 

0.2 

0.2 

0.2 

0.2 
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Table  5.2 


H  =  2  km 

H  =  8 

H  =  32 

CO 

OJ 

1 — 1 

II 

K 

H  =  512 

Y 

Fa 

O' 

0 

0 

0 

0 

0 

1 

3714xxx 

I362XX 

237x 

36 

1 

3 

■■  1095  XXX 

2445 XX 

662x 

109 

2 

10 

ll66xx 

9i6xx 

i494x 

353 

6 

30 

1327X 

1286x 

859x 

842 

18 

1° 

3366 

3333 

2957 

936 

33 

3 

396 

394 

385 

309 

62 

10 

43.5 

43.4 

42.9 

40.5 

28.0 

30 

7.5 

7.5 

7.5 

7.1 

5.9 

6o 

1.2 

1.2 

1.2 

1.2 

1.1 

90 

-  1.4 

-  1.4 

-  1.3 

-  1.3 

-  1.0 

120 

-  2.2 

-  2.2 

-  2.1 

-  2.0 

-  1.7 

150 

-  1.5 

-  1.5 

-  1.5 

-  1.4 

-  1.2 

i8o 

0 

0 

0 

0 

0 

O' 

0 

0 

0 

0 

0 

1 

3713XXX 

1360XX 

229x 

36 

1 

3 

1093XXX 

2442XX 

660x 

107 

2 

10 

ll60xx 

9l4xx 

i487x 

348 

6 

30 

131  Ox 

1271X 

85  Ox 

827, 

16 

1° 

3280 

3252 

2894 

915 

31 

'3 

364 

363 

357 

289 

72 

10 

32.7 

32.7 

32.5 

31.2 

22.3 

30 

3.6 

3.6 

3.6 

3.5 

3.1 

6o 

0.9 

0.9 

0.9 

0.8 

0.8 

90 

0.4 

0.4 

0.4 

0.3 

0.3 

120 

0.3 

0.3 

0.3 

0.3 

0.3 

150 

0.3 

0.3 

0.3 

0.3 

0.2 

l8o 

0.2 

0.2 

0.2 

0.2 

0.2 
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(5.29) 

^  =  (V  ^p^  ' 

(5.30) 

n  =  —  +  —  (m^  +  JL 
R  2 

(5.31) 

A®  =  +  m®  + 

(5.32) 

n 

(5.33) 

m 

fa  cos  Q:  =  ^ 

(5.3J+) 

i 

fa  sin  q:  =  ^ 

These  formulas  can  he  used  only  when  'i'  <  20°.  On  the  other  hand, 

the  effect  of  zones  Y  >  20°  is  often  negligible.  Otherwise,  this  effect 

may  be  computed  for  a  few  points  only  and  interpolated  for  the  other 

points.  When  the  components  5,6,  6. have  been  computed,  we  must 

n  m  t 

convert  them  to  the  system  x,y,z  by  the  matrix  multiplication: 

/  5 

/  ^  1  /  '  ' 

(5.35)  I  5  I  =  I  -  cos  cp  sin  X  -  sin  cp  sin  X  cos  X 

cos  cp  0 

The  components 


/  5  \  / 

X 

5 

— 

y 

\\l 

\ 

m 


6  along  the  trajectory 

5  horizontally  to  the  right  hand  side 
R 

6^  perpendicularly  (not  vertically)  down  . 

can  be  obtained  as  follows.  Take  three  consecutive  points  P^ 

and  P  of  the  trajectory.  Compute: 
i+1 


2k 


X'  =  X  -  X 

i+1  i-1 


(5.36)  y  =  -  y^_^ 


z '  =  z  -  z 

i+1  i-1 


„  (x'xj+y'yj)  zj-z'pi* 


X  = 


r'p^r^ 


(5.37) 


r  'p 

i 


2"  _  ^'^i+y’y-i+^*^i 

r  'p 

i 


(5.38)  p"  =  l/x"2  + 


n"2 


/S 


A 


R 


D 


/.z" 

0 

\  P" 


-  X 


x”z" 


y"  \ 


zX 

-n"  / 


/  " 


n 


m 


Size  of  Blocks.  Around  the  point  P  define  four  regions,  each  bounded 
by  latitudes  cpi,  cpg  and  longitudes  \g  . 


Region  A 


'te-  Ti  =  3 

Xs-  >1  = 

Icp  -  cppi  <  1° 

I  \  -  Ipl  <  1°  30’ 


Region  B , 
outside  A 


cfe-  =  7° 

Xs-  Xi  =  9° 

Icp  -  cppl  <3° 

IX  -  Xpl'  <  4° 
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Region  C; 


cfte-  cPi  =  25® 
outside  B  ^2-  =  30° 

I  CD  -  CD  I  <10 
^  P 

IX  -  X  I  <  12°  30' 

P 

^1}  2-11  divisible  "by  5®. 

Region  D  outside  C 

For  region  A  ,  pick  up  the  mean  anomaly  cards  5'><  5'  (3^  X  48  = 

1728  cards),  for  region  B  ,  20'X  20'  cards  (21  X  27  -  9  X  12  =  459  cards) 
for  region  C  ,  l°x  1°  cards  (25  X  30  -  7  X  9  =  687  cards)  and  for 
region  D  ,  5°x  5°  cards.  Region  D  will  be  used  for  a  few  points 
only. 


The  influence  of  zone  D  on  5  , 

X 


by  means  of  the  Cartesian  coordinates 


x,y,z  are  computed  for  the  centers  of 


integration  can  be  carried  out  by  the 


(5.40)  =  aa 


D 


D  ^ 

O  Z  -  Z 

Ac  R  r  M  p  ^ 

/as  =  —  P  •  — Ta —  da 
z  2tt 

D  ^ 


5,6  can  also  be  computed  directly 

y  z 

x,y,z.  If  these  coordinates 
the  5°X  5®  squares,  then  the 
formulas 


where 


(5.4i)  =  (x  -  X  )^  +  (y  -  y  )‘^  +  (z  -  z  )^ 

^  ^  M  P  M  P  M  P 

and  As  ,  As  ,  As  are  the  contribution  of  the  zone  D  to  5  ,  5  ,  6  . 

X  y  ^  z  X'  y  z 
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PART  II 
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We  shall  assume  uniform  coverage  of  the  whole  earth  by  gravity 


stations  so  that  c  approximately  depends  on  the  relative  position  of 
points  P(YjO:)  and  P'(Y',a')  only.  Fiurthermore  we  assume  that  the  error 
function  has  a  sharp  maximum  at  P'  =  P  and  drops  off  rapidly  to  zero 
with  increasing  distance  PP’ . 

Then  we  can  simplify  the  integrations  considerably  by  approximately 
replacing  a(Y,o:,  Y' ^o:' )  by 

(6.2)  •  6(Y'  -  Y)  6(a'  -  a) 
sin  I 

where  6(Y'  -  Y)^  6(a'  -  a)  are  Dirac's  delta.’ functions  and  S  is  a  constant 
given  by 

2tt  it 

(6.3)  S  =  /  J  CT(Y,a,Y',Q:')  sin  Y'  dY'  da'. 

a'=0  Y'Lo 

Since  according  to  a  property  of  the  delta  function 


J'f(x')  6(x'  -  x)  dx'  =  f(x), 

the  integrations  ■'.■ri.th  respect  to  Y',  a'  can  be  performed  immediately 
and  we  get 


q  2"  " 
0=0  Y=0 


-  t  cos  Y)®  sin  Y  dY  da, 


m 


m 


S 

4tt® 


2tt  n  6  ('cos^a'j 

r  r  — g-  sin^Y )  8  I 

D  I  sin®aj  " 

0=0  Y=0  ^ 
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and  finally, 


(6.4) 


rr,®  — 

m  =  — 
n  2tt 


t*(t^+l)  _jb 
2(t*-l)* 


t®-l 


a  a 

=  ni-  =  TT’ 
m  £  2tt 


t‘^(t^+l)  .  |t- 

4(t*.i)®  8  ^"rt+: 


Since 


4-  ^  1  H 

t  =  -  =  1  _  —  +  . 

R  +  H  R 


(Rrmean  radius  of  the  earth)  we  can  develop  in  series,  and  neglecting 
H 


higher  powers  of  ~  we  get 
R 


(6.5) 


rv,® 

m  = 


SR' 


n  8ttH' 


a  ) 


m®  =  m® 


SR' 


m 


JL  16ttH®  ■ 


(For  H  =  500  km  the  error  of  these  approximations  is  smaller  than  7^.) 
Validity  of  the  Simplified  Integration 

Now  we  have  to  investigate  in  more  detail  the  validity  of  replac¬ 
ing  the  error  function  a  by  a  product  of  Delta  functions.  For 
simplicity  we  limit  ourselves  to  the  plane  approximation.  In  this 
case,  eqs.  (5 .20)  -  (5.22)  are  simplified  to 
H 


6  = 
n  2t! 


•r-T 


(6.6) 


6 

m 


=  ^  J’J'  ' 


where  the  xy- plane  is  horizontal  and  the  z-axis  contains  point  P. 
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These  equations  have  the  general  fo?rm 
00 

(6-7)  ;;  u(x  ^y)  f(x,y)  dx  dy 

•00 

and  the  mean  square  error  is  given  hy 


(6.8) 


m 


=  M  a(x,y,x',y')  f(x,y)  f(x',y')  dx  dy  dx'  dy'  . 


As  before  we  assume  the  error  function  CT(x,y^x',y')  to  have 
appreciable  magnitude  only  for  (x* ,y ' )=(x,y)  and  to  be  almost  zero 
elsewhere.  Then^  only  the  points  (x' ,y ' )=(x,y)  will  contribute  signi¬ 
ficantly  to  the  integral  and  if  we  put 
x'=x+§,  y'=y  +  '[], 

then  only  small  values  T]  need  to  be  considered.  If  the  usual 

conditions  of  continuity  and  differentiability  are  satisfied, 

f(x',y')  can  be  developed  in  a  Taylor  series  with  respect  to  I,  T]: 

f(x',y')  =  f(x,y)  +  §  f  +  71  f  + 

X  y 

+  i  §®f  +  §7]  f  +  I  71=f  +  .  .  . 

2  XX  xy  2  yy 

where  f  ,  f  ,  .  denote  partial  derivatives  at  the  point  (x,y)  and 

X  XX 

5  =  x'-x,  71=y'-y. 

Inserting  this  in  (6.8)  we  can  separate  the  integrations  over  x,  y 
and  over  x',  y',  getting 
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m 


=  JT  ^(x,y,x',y')  dx’  dy'  •  [  f(x,y)]  *  dx  dy  + 

-00  -00 
00  00 

+  JJ  (x'-x)  CT(x,y,x',y')  dx'  dy'  •  dx  dy  + 

-00  -00 

00  00 

+  /J*  (y'-y)  ^(x,y,x',y’)  dx'  dy'  •  JJ  ff  dx  dy  + 


+  p  JJ*  (x'-x)®  a(x,y,x',y')  dx'  dy'  •  JJ  dx  dy  + 


+  (x'-x)(y'-y)  a(x,y,x',y')  dx'  dy'  •  JJ  ff  dx  dy  + 


-00 

00 


xy 


,y')  dx'  dy  '  '  n  ^  ^  i 

M  ■/  1/ 


provided  aix,y,x',y')  depends  on  x-x',  y-y'  only,  i.e.,  we  have  the 
same  accuracy  everywhere.  Furthermore,  if  o  is  a  symmetrical 
function  of  x'-x  and  y'-y  which  is  a  very  natural  assumption,  then 
the  integrals  containing  x'-x  and  y'-y  linearly  will  vanish  and  there 
remains 


m 


=  So  J’J’[f(x,y)]®  dx  dy  + 


(6.9) 


00  '00 

+  I  Si  JJ  ff^  dx  dy  +  I  Ss  J’J*  ff^  dx  dy  . 
—  00  —00 


where 


So  =  JJ  a  dx'  dy'.  Si  =  /]“  (x'-x)®  odx'  dy'. 


(6.10) 


—  CO 
00 


Sz  =  JJ  (y'-y)®^  <ix'  dy' 
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If  we  would  have  replaced^  as  before,  the  error  function  ct  by  a  prod¬ 


uct  of  delta  functions, 

So  •  6(x'-x)  B(y’-y)  , 
we  would  have  got 

00 

JJ  [f(x,y)]®  dx  dy  , 

-00 

which  Is  the  first  and  principal  term  of  (6.9). 

Thus  the  use  of  delta  functions  Is  justified.  Transferring  these 
results  from  the  plane  to  the  spherical  case  offers  no  difficulties. 
Sq,  the  integral  of  the  error  function  over  the  plane  (6.IO),  is 
replaced  by  S,  the  integral  over  the  unit  sphere  (6.3)j  we  have  the 
relation 

So 

(6.11)  S  =  — 

R® 

where  the  scaling  factor  R^  represents  the  transition  from  the 
terrestrial  sphere  (or  its  tangential  plane)  to  the  unit  sphere. 

The  error  of  the  delta  function  method  is  the  effect  of  Si  and 
S2  in  (6.9). 

In  order  to  evaluate  the  quantities  S,  So,  Si,  Sg  we  have  to 
make  some  assumption  concerning  the  error  function  o  .  If  this 
function  has  the  character  assumed  above  -  everywhere  the  same  form, 
sharp  maximum  at  x'  =  x,  y'  =  y,  vanishing  farther  away  from  this 
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point  -  then,  within  the  limit  of  the  obtainable  ac:cura';y ,  v/e  ''a/j 


usually  take  the  well-knovm  function 

(  1  i\  -c*[(x'-x)®  +  (y'-y)®] 

(6.12)  a(x,y,x',y')  = 

(6.10)  can  then  be  evaluated  immediately  to  give 


TTGo 


rra„ 


(6.13)  So  =  -  ,  Si  =  S2  =  - 


2c" 


and 


(6.14)  S  = 


c®R® 

Nov  ve  can  return  to  eqs.  (6.6).  First  we  consider  S  .  Here 

n 


f(x,y)  -  (H®+x®+y®)®^®  * 

Differentiating  twice,  inserting  in  (6.9)  and  integrating  yields 

So  SCSi+Sg) 

=  -  -  -  . 

8ttH*  64 

The  first  terra  is  clearly  the  same  as  the  first  equation  in  (6.5). 
Inserting  (6.13)  we  finally  get 

(6.15)  ra®  =  -  (1  -  — ^  +  .  .  .  ) 

8cV  8c^H® 


and  in  a  similar  way  we  find 


(6.16)  m®  =  =  -  (1 

”  *  i6==ir' 


3 

8c®h® 


.  .  ) 


The  correction  term  in  the  parenthesis  is  practically  zero  for 
elevations  H  »  c  ^  ,  but  in  view  of  the  practical  purpose  of  these 


33 


-1 

accuracy  formulas,  it  can  already  be  neglected  for  H  =  c 
Discussion 

To  sum  up  the  formulas  we  have  for  the  standeird  errors  of  S  , 

n 


(6.17) 


m  =  m, 
m  i 


kg 

H 


where 


(6.18)  kj 


These  formulas  are  valid  for  elevations  H  >  c  .  For  smaller 
elevations  they  cannot  be  applied  (for  H  =  0  they  would  yield  *  , 
which  is  obviously  wrong). 

The  error  covariance  function  o  is  supposed  to  have  the  form 


(6.12)  ; 

-c®s^ 

(6.19)  a(x,y,x’,y')  =  Ooe 

where  s  is  the  distance  of  the  points  (x,y)  and  (x',y')»  According 

to  the  derivation  of  the  above  formulas,  a  is  the  error  function  of 

37 

p,  ,  the  density  of  the  coating.  Since  p=-Ag  +  TT  N  ,  it  contains  the 

2R 

influence  of  errors  in  Ag  and  in  N  .  The  errors  of  N  have  a 
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small  influence  but  they  diminish  much  less  with  increasing  distance 
and  can  hardly  be  represented  in  the  form  (6.19).  So,  in  order  to 
apply  our  formulas,  we  must  be  able  to  neglect  the  errors  of  N.  That 
this  is  possible  can  be  seen  in  the  following  way. 

We  would  have  avoided  the  trouble  concerning  the  errors  of  N  , 
had  we  used  the  direct  formulas  (5.15)  -  (5.19)  instead  of  using  the 
coating  method  (5.20)  -  (5.214-).  These  formulas  are,  however,  very 
difficiilt  to  handle  for  our  purpose.  On  the  other  hand,  for  smaller 
values  of  ^  ,  where  the  effect  of  errors  in  [j,  is  largest,  F  and 
f  are  closely  equal  so  that,  for  this  particular  purpose,  we  might 
replace  p,  in  (5.20)  -  (5.214-)  by  Ag.  Then  we  can  take  a  in  (6.I) 
and  in  the  subsequent  developments  to  be  the  error  covariance  function 
of  the  gravity  anomalies  Ag  only.  From  this  we  can  conclude  that 
the  error  of  (6.17)  due  to  neglecting  the  inaccuracy  of  W  must  be 
small. 

Thus,  we  can  consider  <7  in  (6.19)  to  be  the  error  function  for 
the  interpolation  of  gra’/ity  anomalies.  Then, 

3 

Co  — 

is  the  square  of  the  average  standard  error  of  interpolation,  m  . 

In  Table  6.1  some  nimierical  values  of  m  and  c  are  given,  together 
with  the  constants  k^  and  kg  defined  above. 
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Table  6.1 


Constants  k^ , 

k  2  for  average 

station 

distances  s 

s 

m 

-1 

c 

kg 

km 

mgal 

km 

mgal .km 

mgal. km 

10 

1 

13 

2.5 

1.8 

50 

9 

25 

45 

32 

100 

Ik 

45 

126 

89 

200 

16 

90 

287 

203 

To  get  the  standard  errors  of  6'  and  of  6'  ,  6  ^  the  above 

n'  m  £ 


values  for  and  kg  have  to  be  divided  by  the  elevation 
H  in  kilometers. 

For  the  computation  of  the  dependence  of  m  and  c  on  the 
average  station  distances,  an  idealized  gravity  station  net  consisting 
of  equilateral  triangles  with  side  s  has  been  assumed.  The  com¬ 
putation  was  done  in  a  way  described  in  [Moritz ,  1962b,  Appendix], 
using  numerical  values  for  the  covariance  function  of  gravity  anomalies  given 
in  [  Kaula,  1959]  combined  with  the  standard  interpolation  errors  of 
[Molodenskii  et.  al.,  19^2,  p.  172].  Since  the  method  used  is  not 
quite  rigorous,  the  resxilts  of  Table  2  must  be  considered  somewhat 
preliminary. 
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Error  Correlation  of  6  ,  5  ,  5 

n'  m  SL 


For  the  error  covariances  of  &  ,  S  ,  S .  in  two  points  P,  P' 
-  n  m 

we  have  found  [Moritz,  1963]  : 


n,PP' 


°  n  3  ^  1 

U  -  8  jj8  J  ^ 


8c®h® 


(6.20)  a  '  = 


16^?  '"  ■  3  F  ^  *  .  '  ’ 


i,pp' 


3  s® 

-  [1  -  i®  ^  Sin®a)  +  .  .  .  ]  . 


16c®H® 


P,  P*  are  assumed  to  have  the  same  elevation  H,  s  is  the  distance 
and  OC  is  the  azimuth  fran  P  to  P' .  These  formulas,  being  the  first 
terms  of  power  series,  are  valid  for  small  distances  s  only.  For 
P  =  P*,  by  (6.5)  and  (6.l4)  , 


(6.21)  a  =  ra®,  a  =  m®,  a  =  m®  : 

''  '  '  n,PP  n'  m,PP  m'  £,PP  £,  ' 

the  error  covariances  become  variances,  i.e.,  the  squares  of  the 

standard  errors  of  6  ,  5  ,  5.  , 

n  m  £ 

Dividing  the  covariances  (6.20)  by  the  variances  (6.21),  we  get 


the  correlation  coefficients 


P 


(6.22) 


n,PP’ 


^m,PP' 


ji,PP' 


1  3  s® 

=  1-0  TTa"  + 


=  1  - 


=  1  - 


8  H® 

3  s® 

—7  — 5-  (1  +  2  cos®a)  +  .  . 
16  H®  ^  ‘ 


16  ^  ^  Ein®0!)  +  .  . 
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p  =  Jb  1  means  maximum  correlation,  i.e.,  complete  functional 
dependence;  p  =  0  means  independence  (more  strictly  speaking,  com¬ 
plete  lack  of  correlation) .  For  a  definite  s  ,  the  correlation 
coefficients  (6.22)  will  be  the  closer  to  1,  the  greater  the  altitude 
H  is.  Thus,  with  increasing  altitude,  the  standard  errors  of  6^, 

5,5.  decrease,  but  the  correlation  increases.  Since  the  accuracy 
m  H  - 

is  characterized  by  the  standard  errors  and  the  error  correlation, 
this  is  significant,  too. 

7.  INFLUENCE  OF  THE  REMOTE  ZONES 
The  purpose  of  the  estimation  of  the  influence  of  the  distant 
zones  is  to  decide  how  far  the  integration  must  be  extended  if  certain 
accuracy  requirements  are  prescribed. 

Since  the  effect  of  the  distant  zones  is  almost  independent  of  the 
elevation  H  ,  we  may,  for  simplicity,  put  H=0,  i.e.,  we  consider  P 
to  be  situated  directly  on  the  geoid.  Then,  since 

8^  =  Ag  +  K,  =  r?,  h  =  - 

we  can,  for  the  direct  method,  use  one  of  the  evaluations  of  the  effect 
of  distant  zones  on  the  conponents  §  ,  T]  of  the  deflection  of  the 
vertical  and  on  N  (e.g.,  [Raula,  1957])*  For  our  purpose  it  is  best 
to  use  the  formulas  of  [Molodenskii  et  al.,  I962]  which  can  also  be 
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easily  adapted  to  the  cpating  method. 

Direct  Method  (using  Ag).  The  disturbing  potential  is 


2n  n 

(7.1. a)  T  =  ‘yN  =  —  J  J  ^g  S(cos  Y)  sin  7  d7  da 

^  0=0  ^=0 


where  S(cos  Y)  is  Stokes’  function 

1  Y  Y 

(7. 2. a)  S(cos  Y)  =  - y  -  3  cos  Y  jJn  (sin  —  +  sin®-^)  - 

sin  - 
,  Y 

-  6  sin  —  +1-5  CO  sY. 


If  we  extend  the  integration  .with  respect  to  Y  only  up  to  a  spherical 

radius  Y  <  rr,  we  commit  an  error. 

0 

2n  n 

(7. 3. a)  An  =  —  I  J  Ag  S(cos  Y)  sin  Y  dY  da 
^0=0  Y=Y 

o 

which  can  he  also  written 


„  2tt  n 

(7. a)  An  =  r — •  J*  J  Ag  S(cos  Y)  sin  Y  d'-  da 
^0=0  Y=o 


where 

(7.5 .a)  S(cos  Y) 


rO  if  Y  <  Y 

I 

L  S(cos  Y)  if  Y  >  Y 


Expand  S(cos  Y)  in  a  series  of  Legendre’s  polynomials: 
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03 


(7. 6. a)  S(cos  Y)  =  r 


n=0 


— —  Q  P  (cos  'K) 
2  n  n 


where  the  coefficients  Q  depend  on  Y  .  Then,  by  multiplying  both 

n  o 

sides  by  P  (cos  Y)  sin  Y  and  integrating  from  0  to  rr  we  find 


TT 


TT 


Q  =  r  S(cos  Y)  P  (cos  Y)  sin  Y  dY  =  f  S(cos  Y)  P  (cos  Y)sin  Y  dY  . 

n  w  n  n 


substituting 


(7.7-a)  z  =  sin  —  ,  t  =  sin  — “ 


I. 

2 


we  get 


(7. 8. a)  Q  =  -4  f  P  (1  -  2z®)  S(l  -  2z^)  z  dz  . 
n  .  n 


Performing  the  integration  yields  (we  need  the  only  for  n  ^  2): 

=  2  -  4t  +  5t^  +  l4t^  -  ^  t^  .  30t^  +  47t^  +  l8t^  _  ^  t®  + 

+  (6t^  -  24t^  +  36t^  -  l8t^)  in  t(i+t)^ 


(7. 9. a)  Q, 


^  -  1  -  4t  +  5t^  +  22t^  -  46t^  -  ^  t^  +  136t*^  -i-  1041"^  - 

-  i66t^  -  48t^  +  t^^  +  (6t^  -  42t^  +  io8t^  -  I20t^  + 

5 

+  48t^^)  in  t(l+t)  , 


2  2  Q8  ^  *5  8  T 

=  -  -  4t  +  5t  +  ^t^  -  72t'  -  I56t^  +  320t  ■+  36ot'  _645t 

1120^9  .10  ,11  12  ,,  2  .,4 

- - 1  +  6o2t  +  i4ot  -  2iot  +  (6t  -  66t  + 

3 

8  8  10  ~|  2 

+  26ot  -  48ot  +  420t  -  l40t  )  in  t(l+t). 


8 
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Higher  order  Q  (up  to  8th  order)  may  he  found  in  [Molodenskii  et  al.,  1952 , 
n  -  ’ 

p.  149]. 

Hy  (T*4.a)  and  (7. 6. a)  we  get 
00 

(T.10.'.a)  An  =  ^  S  Q  Ag 

27  „  n  n 


if 

n 

Ag  =  s  Ag 

n 


is  the  development  of  Ag  in  Laplace's  spherical  harmonics  (  As 

o 

Agi  8ire  missing,  as  usual). 

The  coating  method  can  be  treated  in  an  exactly  analogous  way. 
Coating  Method  (using  tJi).  The  disturbing  potential  is 

2n  TT 

(T.l.b)  T  =  'yN  =  'j^  J  J  M«M  (cos  Y)  sin  Y  dY  dCC 

0=0  Y=0 

where 

(T.2.b)  M(cos  Y)  =  * 

sin  - 

If  we  extend  the  integration  with  respect  to  Y  only  up  to  a  spherical 

radius  Y  <  rr,  we  commit  an  error 
o 

2T7  tt 

(T.S.'b)  An  =  J  J  M-  M  (cos  7)  sin  ¥  dY  do 
^0=0  Y=Y 

o 


4i 


which  can  he  also  written 


R  " 

(7.4,h)  An  =1^  J  J  (cos  Y)  sin  Y  dY  da 


2tt  n 


■  rry 


0=0  Y=0 


where 


(7.5.h)  M(cos  y)  =J 


ro 

M  (cos  Y) 


sin 


if  Y  <  Y 


if  Y  >Y 


Expand  M(cos  Y)  in  a  series  of  Legendre’s  polynomials: 

CO 

(7.6.h)  M(cos  Y)  =  Z  (n-l)  q  P  (cos  Y) 

^  n  n 

n=:0 

where  the  coefficients  q  depend  on  Y  .  Then,  hy  multiplying  hoth  side 

n  o 

hy  P  (cos  Y)  sin  y  and  integrating  from  0  to  tt  we  find 
n 


q  = 
n 


:  7 - .  r  M  (cos  Y)  P  (cos  Y)  sin  Y  dY  = 

2(n-l)  J  '  '  n  ' 

=21^1)  5  ^  "n  ''  ""  • 

Y  sin— 
o  2 


;5y  substituting 

Y  Y 

(7. 7. a)  z  =  sin  t  =  sin 

we  get 

(7.8.h)  q  =  J  P  (l-2z®)  dz. 

n  n-l  'J  n 
t 
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Performing  the  integration  yields  (we  need  the  q  only  for  n  ^  2): 

n  — 

q^  =  |(l  -  7t  +  28t^-  42t^+  20t^), 

q^  =  |(1  -  9t  +  60t^-  l62t^+  l80t^-  70t^), 

q^  =  •|(1  -  lit  +  liot^-  462t^+  880t'^-  770t^+  252t^^), 

(7,9.h)  q^  =  |(1  -  I3t  +  l82t^-  1092t^+  3l20t^-  4550t^+  3276t^^-  924 t^^), 

q^  =  |(l  -  15t  +  280t^-  2268 t^+  9000t^-  1925 0t^+  22680t^^- 

-  I3860t^^+  3432t^^), 

q-  =  -(1  -  17t  +  4o8t^-  4284 t^+  2244ot^-  6545 0t^+  lll384t^^- 
8  7 

-  109956t^^+  58344 1^^-  128701^’^)  . 


By  (7.4.h)  and  (7.6.8)  we  get 


*  R  n-1 

An  =  -  S 


y  2n+l  ^n^n 
n=2 


if 


n 

IJi  =  I!  H. 


n=2 


n 


if  the  development  of  p  in  Laplace's  spherical  harmonics  (p  and 

o 

are  missing  together  with  Ag  Agi). 

o 

If 


Ag  =  2  Ag 


n=2 


n 
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then 


.  =  2  ,  As. 

y  ^  n-l 

n=2 


Therefore, 


11  = 


^  ‘  n-2 


SO  that 


2n+l 


n 


-|  \ 

2(n-lj  n 


Inserting  this  in  the  formula  for  N  we  finally  obtain 


E 


(T.io.h)  An  =  —  Z  Ag^, 

^  n=2 

which  is  analogous  to  (Y.lO.a). 


Mean  Square  Effect  on  N,  f,  11,  Squaring  (7. 10. a)  and  averaging  over 
the  whole  earth  we  get 


_ _ 

AN=  =  — 2  E  S  Q  Q  ,  •  Ag  Ag  , 

ky  ,  „  n  n'  n  n' 

n=2  n'=2 


Since  the  integral,  over  the  sphere,  of  the  product  of  two  Laplace 
harmonics  of  different  degree  is  zero  (Orthogonality  l),  only  the 
terms  where  n'  =  n  remain  and  we  obtain 


TT-a  R' 


(7. 11. a)  An®  = 


4y®  ^ 

n=2 


*  Ag  ® 

n  n 


where  AN  is  the  root  mean  square  influence  on  N  of  the  zones  beyond 


a  spherical  radius  Yq  ,  and Ag  is  the  root  mean  square  average  of 
the  Laplace  spherical  harmonic  Ag 
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The  derivation  of  the  RMS  influence  on  the  components  Tl  of  the 
deflection  of  the  vertical  is  rather  involved.  Since  it  is  omitted  in 
[Molodenskii  et.  al.,  I962],  it  may  he  of  interest  to  give  it  here. 


By  (T.lO.a), 


SAg 


(7.12) 


Ati  = 


R  a«P  2y  n  dcp 


1  3  K 


3Ag 


S  Q 


R  cos  cp  3\  2v  cos  cp  ^  n  31 
^  ^  n=2 


for  the  direct  method.  (For  the  coating  method  replace  Q  hy  q  . ) 

n  n 

Hence, 


^  gn  S  gn' 


^  ^n\'  3cp  3cp  ^ 

^  n=2  n'=2  't' 


2 _  5  g^  3  gn' 


1  2  ^  ^  Q  Q  ,  s  V,  f 

4V^  „  ,  ^  n  n'  cos  cp  31  31 


n=2  n'=2 


so  that 


GO  00 


5  gr.  5  g^,  1  3  3  gn. 


(7.13)  +  AT1»  =  E  E  Q  Q  , 

4'y  n  n  3cp  3cp  Cos  cp  31  31 

^  n=2  n'=2  'r 

(The  bar  denotes  averaging  over  the  sphere.) 


Ag  can  be  written  as 


m 


(7.14)  Ag  =  2  (a  cos  ml  +  b  sin  m\)  P  (sin  cp) 

n  nm  nm  n 

m=0 


where  a  ,  b  are  constant  coefficients  and  cp,  1  are  geographical 
nm  nm 


coordinates.  Now  we  can  differentiate  with  respect  to  co  and  1: 


(7.15) 


Hence  we 


(7.16) 


Here  we 


(7.17) 


m 


c|Ag  n  dP 

n  _  ,  ,  ,  ,  n 

— r —  =  Z  (a  cos  niA.  +  b  sin  mAj  -  , 

ocp  nm  nm  dcp 

m=0 


SAgn  n  ^ 

— r- —  =  E  (mb  cos  mA  -  ma  sin  mA)  P 
oA  nm  nm  n 

m=0 


get,  taking  n'  ^  n  , 


4rr 


^6  ,  2tt  tt/P  5Ag  SAg 


n  n 


;  j 


n  n' 


‘  aIo  cp=  -V2  ^'P 


cos  cp  dcp  dA 


=  2tt  a  a  .  f 

n  *  r\  w 


tt/2  dP  °  dP  ° 
n  n' 


no  n'o  dcp  d9 


cos  cp  dcp  + 


n  tt/2  m  m 

^  ,  N  P  dP  dP  , 

+  TTZ(aa  +bb,)  n  n'  cos  cp  dcD 

^  nm  n’m  nm  n'm  .  — - — 

m=l  -n/2  dcp  dcp 


have  used  the  well-known  orthogonality  relations 


2tt 


J  cos  pA  cos  qA  dA  =-| 
0 


r  0.  p  q 

rr,  p  =  q  0  , 

2Tr  ,  p  =  q  =  0 


2rr 

; 


cos  pA  sin  qA  dA  =  0  , 


2tt  7  0,  p  7^  q 

J  sin  pA  sin  qA  dA  =  <  n,  p  =  q  7^  0 
0  lo,  p=q=0 
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Similarly 


where 


m  ,  V  (n+iri-2)l  ,  (n+ni-4)l 

C  =  2(2n  -  1)  +  2(2n-5)  j ^  +  .  .  .  , 

n  (n-m.)!  (n-ni-2)L 

,  m  ^  It  i^ti 

(Concerning  C  ,  there  is  a  misprint  on  p.  165,  op.  cit.;  +  C 
n  n 

should  read  =  C  ' . ) 

n 


Hence  we  findj  for  n'  ^  n: 


n 


m. 


♦■TT 


aAg  SAg 

_ n  _ n' 

Sep  9cp 


n  E  (a  a,  +b  h)»(  -m*C  ) 
nm  n'm  nm  n'm  n  , 

Itb=l 

n  -n  even^ 


0^  n'  -  n  odd, 


n 


m 


4tt 


1 


SAg  3Ag 


n  n' 


cos^cp 


^cp  9X  SX 


n  S  (a  a  .  +  b  b  ,  )  ‘m^C  , 

nm  n'm  nm  n'm  n 

m=:l 

n'-n  even, 


,  0,  n’  -  n  odd 


so  that 


2lAg  9Ag  9^  9Ag  , 

n  n'  1  n  n' 

(T«19)  ^  '  - 2  TT 

ocp  ocp  cos  cp  oX  oX 
for  n'  7^  n. 


r-  n 


For  n'  =  n  we  find  in  the  same  way 


/SAg  \ 

n 

2 

2 

■  n(n+l) 

9cp  i 

2  1 
COS  cp  ’ 

ax  J 

2n+l 

3  ^  1  niwli  I  l£l4l(3  =  +  b  p 

no  2  2n+l  m-1  (n-m)i  nm  nm 
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Forming  the  average  of  Ag  ®  itself  we  find  easily 

n 

_ 1  2Tr  Tr/2 

=  W  J*  J*  tP  'i'P  = 

\=0  cp=  -tt/2 


(7.20) 


,  ,  ,  ,  ,  r  (a  >  *  t  =) 

2n+l  no  2(2n+l)  .  (n-m)  I  run  nm 

m=l 


2  1 

a  + 


so  that 


(7.21) 


2 

1 

\ 

+  a 

cos‘’cp 

} 

=  n(n+l)2lg 


Taking  (7. 19)  and  (7.21)  into  account,  we  find  hy  (7.13) 

00 

A?®  +  =  1  ^  n{n+l)  Q  . 

4v  n  n 

'  n=2 

For  reasons  of  symmetry,  and  All  can  be  taken  to  be  equal,  and  we 
finally  get 


(7. 22. a)  =  Si®  =  E  n(n+l)  Q_®  Ag_® 


■n-O 


n  n 


for  the  direct  method. 


For  the  coating  method  we  must  replace  Q  by  q  ,  obtaining 
-  n  n 


(7.11.b)  S*  =  ^  Z  q  ®  Ag  , 

4y^  n  n 

'  n=2 


(7.22.b)  Af  =  ATI^  =  ^  2  n(n+l)  q  ®  Ag  ®  . 

oy  n  n 

'  n=2 
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Relationship  ‘between  Ag  and  the  Covariance  Function.  Comparing 
eq.  (7.20)  with  eq.  (33)  of  [Kaula,  19591  and  taking  the  relationship 
"between  conventional  and  fully  normalized  spherical  harmonics  into 
account,  we  find  that 

(7.23)  Ag  =  =  c  , 

n  n 

the  coefficient  of  n-th  degree  in  the  development  of  the  covariance 
function  of  gravity  anomalies,  C(s),  in  a  series  of  Legendre's  poly¬ 
nomials; 

00 

(7.24)  c(s)  =  E  c  P  (cos  s). 

„  n  n 


Mean  Square  Effect  on  6  ,  6  ,  5  .  In  order  to  avoid  confusion,  we 

- K  —  M  —  L 

shall  change  our  notation  of  the  components  of  the  gravity  disturbance 
vector,  now  writing  as  subscripts,  capital  letters  N,  M,  L  instead 
of  n,  m,  Jl. 

Then, 

(7.25a)  5^  =  Ag  +  ^  N,  6^  =  y|,  5^  =  yT]. 

Denoting  the  RMS  effect  of  the  zones  beyond  7  on  these  components 

o 

by  As  ,  As  ,  As  ,  respectively,  we  find  from  (7.11.a,b)  and  (7.22.a,b), 
IT  M  L 

taking  (7.23)  into  account,  the  following  formulas; 
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Direct  Method 


55:*  =  2  Q_ 


N 


n=2 


n  n 


(7. 26. a) 


Ad  =  Ae”®  =  o  2  n(n+l)  Q  ®  c 
M  L  o  ^  n  n 

n=2 


If  we  use  the  coating  li  ,  then 


(7.25.b)  S^.  ^  + 8^=y?,  =  yH, 


and  we  have  instead: 


Coating  Method 


Aa  ®  ~  E  q  2  c  , 

N  16  ^  n  n 

n=2 

(T.26.h)  « 

AHT,®  =  Q  r  n(n+l)  q  ®  c  . 


M 


8 


n=2 


n  n 


Numerical  Results.  First  we  give  tables  for  the  C 


taken  from  [Molodenskii  et.al.,  1962])^  and  for  q 


computed  by  eqs.  (T-9*b)). 


(Table  7.1»aj 
(Table  7.1.b, 
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Table  7.1. a 


Coefficients  Q  for  the  direct  method  as  functions  of  radius  7 
n  o 


o 

Q 

^2 

s 

% 

% 

0° 

+2.000 

+1.000 

+0.667 

+0.500 

+0.400 

+0.333 

+0.286 

5° 

+1.801 

+0.802 

+0.470 

+0.304 

+0.206 

+0.141 

+0.095 

8.99° ( 1000km) 

+1.634 

+0.639 

+0.312 

+0.152 

+0.061 

+0.004 

-0.032 

10° 

+1.593 

+0.596 

+0.274 

+0.118 

+0.030 

-0.024 

-0.056 

13.49° (1500km) 

+1.457 

+0.472 

+0.159 

+0.015 

-0.058 

-0.095 

-0.110 

17. 99° (2000km) 

+1.306 

+0.339 

+0.049 

-0.070 

-0.116 

-0.126 

-0.118 

20° 

+1.247 

+0.290 

+0.011 

-0.094 

-0.127 

-0.125 

-0.105 

25° 

+1.133 

+0.204 

-0.044 

-0.116 

-0.120 

-0.093 

-0.057 

30° 

+1.062 

+0.161 

-0.058 

-0.103 

-0.087 

-0.050 

-0.013 

40° 

+1.021 

+0.144 

-0.051 

-0.080 

-0.057 

-0.024 

+0.002 

50° 

+1.044 

+0.122 

-0.094 

-0.116 

-0.067 

-0.006 

+0.033 

60° 

+l.o4o 

+0.024 

-0.188 

-0.134 

-0.009 

+0.069 

+0.063 

70° 

+0.958 

-0.130 

-0.240 

-0.053 

+0.095 

+0.081 

-0.015 

80° 

+0.807 

-0.260 

-0.187 

+0.073 

+0.111 

-0.016 

-0.075 

90° 

+0.636 

-0.306 

-0.068 

+0.128 

+0.023 

-0.076 

-0.011 

100° 

+0.504 

-0.274 

+0.024 

+0.090 

-0.046 

-0.034 

+0.041 

110° 

+0.441 

-0.222 

+0.049 

+0.038 

-0.043 

+0.008 

+0.020 

120° 

+0.431 

-0.208 

+0.047 

+0.028 

-0.035 

+0.010 

+0.010 

130° 

+0.430 

-0.237 

+0.078 

+0.019 

-0.050 

+0.034 

-0.003 

140° 

+0.389 

-0.263 

+o.i4o 

-0.039 

-0.025 

+0.050 

-0.045 

150° 

+0.291 

-0.234 

+0.171 

-0.107 

+0.050 

-0.006 

-0.023 

0 

0 

+0.159 

-0.145 

+0.127 

-0.107 

+0.086 

-0.064 

+0.044 

0 

0 

1 — 1 

+0.045 

-0.044 

+0.043 

-o.o4i 

+0.039 

-0.037 

+0.035 

0 

0 

CO 

1 — 1 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 
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Table  Y.l.b 


Coefficients  q  for  the  coating  method  as  functions  of  radius  Y 
n  o 


Y 

o 

^12 

% 

% 

0° 

+2.000 

+1.000 

+0.667 

+0.500 

+0.400 

+0.333 

+0.286 

5° 

+1.565 

+0.697 

+o.4o8 

+0.265 

+0.179 

+0.123 

+0.083 

8.99° ( 1000km) 

+1.226 

+0.465 

+0.215 

+0.095 

+0.026 

-0.016 

-0.042 

10° 

+1.142 

+o.4o8 

+0.170 

+0.056 

-0.007 

-0.044 

-0.066 

13.49° (1500km) 

+0.858 

+0.222 

+0.024 

-0.062 

-0.102 

-0.119 

-0.121 

17. 99° (2000km) 

+0.512 

+0.009 

-0.128 

-0.170 

-0.173 

-0.156 

-0.129 

20° 

+0.366 

-0.075 

-0.182 

-0.202 

-0.185 

-0.152 

-0.113 

ro 

0 

+0.033 

-0.251 

-0.275 

-0.233 

-0.169 

-0.101 

-o.o4o 

30° 

-0.255 

-0.374 

-0.309 

-0.206 

-0.103 

-0.018 

+o.o4o 

0 

0 

-0.676 

-0.460 

-0.228 

-0.045 

+0.066 

+0.103 

+0.085 

0 

0 

-0.878 

-0.363 

-0.037 

+0.116 

+0.127 

+0.057 

-0.021 

60° 

-0.875 

-0.156 

+0.138 

+0.153 

+0.o4o 

-0.058 

-0.071 

0 

0 

-0.707 

+0.070 

+0.206 

+0.061 

-0.074 

-0.073 

+0.007 

0 

0 

CO 

-0.433 

+0.235 

+0.149 

-0.064 

-0.090 

+0.014 

+0.060 

0 

0 

o^ 

-0.121 

+0.293 

+0.018 

-0.119 

-0.007 

+0.068 

+0.003 

100° 

+0.165 

+0.241 

-0.101 

-0.070 

+0.070 

+0.020 

-0.050 

110° 

+0.376 

+0.115 

-0.147 

+0.026 

+0.060 

-0.048 

-0.011 

120° 

+0.484 

-0.028 

-0.105 

+0.086 

-0.012 

-0.039 

+o.o4o 

130° 

+0.483 

-0.137 

-0.014 

+0.068 

-0.058 

+0.020 

+0.015 

l4o° 

+O.406 

-0.178 

+0.065 

+0.002 

-0.035 

+o.o4i 

-0.029 

150° 

+0.275 

-0.153 

+0.092 

-0.050 

+0.018 

+0.004 

-0.016 

160° 

+0.138 

-0.088 

+0.066 

-0.051 

+0.038 

-0.027 

+0.018 

0 

0 

+0.037 

-0.025 

+0.021 

-0.019 

+0.017 

-0.015 

+0.014 

!-■ 

00 

0 

0 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 
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Now  we  are  going  to  evaluate  formulas  (7. 26. a,!)  for  the  direct 
and  the  coating  method,  for  different  sets  of  c  .  First  we  use  the 
maximum  estimates  for  c  of  [Kaula,  1959^  P»  2^19]: 

c^  =  15,  =  43,  Cj^  =  30,  c^  =  Cg  =  c^  =  Cg  =  25  mgal! 

The  results  are  given  in  Table  7.2.  The  summation  has  first  been 
extended  to  n  =  8;  higher  order  c  have  been  neglected.  In  order  to 
see  the  convergence  of  the  series  (7.26.a,b)  we  have  also  performed 
the  summation  up  to  n  =  5  only. 


Table  7-2, 

RMS  influence  of  the  zones  beyond  a  radius  Y  on  6  ,  6  ,  6  ,  based 

o  N  M  L 

on  the  maximum  estimates  for  the  degree  variances  of  [Kaula,  1959] 


Summation  up  to  8th  order  Summation  up  to  5th  order 


Direct  Coating 


5° 

8.99° (1000km) 
10° 

13.49° (15 00km) 
17.99° (2000km) 
20° 

25° 

30° 

40° 

50° 

60° 

70° 

80° 

90° 

100° 

110° 

120° 

130° 

l40° 

150^ 

160° 

170° 

180° 


Direct 

mgal 


9.4 

10.6 

7.8 

8.1 

7.5 

7.7 

6 .6 

6.7 

5.7 

5.9 

5.4 

5.6 

4.7 

4.8 

4.4 

4.1 

4.2 

3.8 

4.2 

4.0 

4.3 

4.3 

4.1 

4.3 

3.8 

4.2 

3.3 

3.7 

2.7 

3.0 

2.3 

2.4 

2.2 

2.3 

2.4 

2.6 

2.4 

3.0 

2.2 

2.8 

1.6 

2.4 

0.6 

1.1 

0.0 

0.0 

Coating 

M  As  =  A5 

K  M  L 
mgal 


2.0 

9.3 

1.4 

6.0 

1.3 

5.4 

1.0 

4.3 

0.7 

4.3 

0.6 

4.4 

0.7 

4.5 

0.8 

4.7 

1.1 

5.2 

1.1 

4.6 

0.9 

4.0 

0,8 

3.3 

0.6 

3.1 

0.5 

2.8 

0.5 

2.6 

0.5 

2,2 

0.5 

2.2 

0.5 

2.2 

0.5 

2.2 

0.4 

1.8 

0.2 

1.3 

0.1 

0.5 

0.0 

0.0 

As,  =  As^ 

N  M  L 
mgal 


9.3 

10.1 

7.8 

8.1 

7.5 

7.6 

6.5 

6.3 

5.6 

5.2 

5.2 

4.9 

4.7 

4.3 

4.3 

4.0 

4.1 

3.7 

4.2 

3.9 

4.2 

4.1 

4.0 

4.0 

3.7 

3.8 

3.3 

3.5 

2.7 

2.9 

2.3 

2.4 

2.2 

2.3 

2.3 

2.5 

2.4 

2.8 

2.2 

2.8 

1.4 

2.0 

0.4 

0.7 

0.0 

0.0 

As  As  =  As 

N  M  L 
mgal 


2.0 

CO 

CO 

1.4 

5.9 

1.3 

5.3 

0.9 

3.4 

0.6 

2.6 

0.5 

2.8 

0.6 

3.8 

0.8 

4.6 

1.0 

4.8 

1.0 

4.3 

0.9 

3.7 

0.8 

3.1 

0.6 

2.8 

0.5 

2.6 

0.5 

2.3 

0.5 

2.0 

0.5 

2.0 

0.5 

2.1 

0.5 

2.0 

0.4 

1.8 

0.2 

1.1 

0.1 

0.4 

0.0 

0.0 

Fig.  7*1  shows  a  graphical  representation  of  the  first  part  of 


Table  7.2. 


55 


RADIUS  OF  NEGLECTED  ZONE 


One  would  expect  a  rather  monotonous  decrease  of  the  cirrves  from 


0°  to  180°.  The  non-monotonousness  is  especially  apparent  in  the 

coating  method.  It  is  probably  -due  mainly  to  neglecting  terms  of 

higher  than  8th  order;  it  is  still  stronger  if  we  use  only  terms  up 

to  5th  order  as  Table  7*2.  shows.  It  seems,  however,  to  be  hardly 

worth-while  to  go  up  higher  than  to  8th  order,  in  view  of  the  increasing 

complexity  of  the  formulas  for  Q  and  q  and  the  uncertainty  of  the 

n  n 

numerical  material  for  c  .  Since  the  higher  order  terms  can  be  ex- 

n 

pected  to  influence  mainly  the  dips,  it  is  justifiable  to  "bridge" 

them  empirically  in  order  to  get  a  smooth  monotonous  function. 

More  recent  values  for  the  c  can  be  obtained  from  [Kaula,  I961] 

n  - 

and  from  [Uotila,  I962],  Both  give  a  development  of  the  gravity  anomaly 
field  in  fully  normalized  spherical  harmonics.  Writing  (7.20)  in  fully 
normalized  harmonic  coefficients  we  get  [Kaula,  1959^  eq.  (33)] 

n 


E 

m=0 


(a  ®  +  b  ^) . 
nm  nm 


By  this  equation,  the  c  can  be  computed. 

n 

From  [Uotila,  1962,  Table  5]  we  get  (flattening  l/297): 


c^  =  22,  c^  =  42,  Cj^  =  32. 


Only  values  up  to  4th  order  are  given.  The  almost  perfect  agreement 
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with  the  c  of  [Kaula,  1959]  is  accidental. 

n  - 

From  [Kaula,  1961]  we  get  for  the  flattening  1/298.24: 

c^=  O.Y,  *^4~  *^5~  ^7~  *^8~ 

if  we  use  the  flattening  l/29T>  then  Cg  is  changed  to  ag=  I7.  The 
influence  of  the  distant  zones,  according  to  these  values,  is  given  in 
Table  7*3  and  in  Fig.  7*2.  We  see  that  the  effect  of  a  wrong  flattening 
is  much  stronger  in  the  direct  than  in  the  coating  method.  This  can  be 
expected  since,  as  Tables  7*l*a  and  7*l*b  show,  qg  is  smaller  and  changes 
sign  more  often  than  Qg , 
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Table  7.3 


EMS  influence  of  the  zones  beyond  a  radius  on  6  ,  6  ,  6  , 

o  N'  M'  L" 

for  two  different  flattenings,  based  on  a  spherical  harmonics 
development  obtained  by  a  combined  adjustment  of  gravimetric,  astro- 
geodetic,  and  satellite  data  [Kaula,  I961]  (summation  up  to  8th 
order) 


flattening  I/298.2U  flattening  1/297 


Direct 

Coating 

Direct 

Coating 

SC 

=  5C 

=  Ab 

=  Ab 

0 

N  M 

L 

N  M 

L 

N  M 

L 

N  M 

L 

mgal 

mgal 

mgal 

mgal 

5° 

h.3 

5.7 

0.9 

5.0 

8.5 

8.6 

1.9 

7.4 

8.99°(lOOOkm) 

3.3 

4.1 

0.6 

3.0 

7.5 

7.1 

1.4 

5.3 

10° 

3.1 

3.8 

0.5 

2.6 

7.2 

6.8 

1.3 

4.8 

13.  (1500km) 

2.5 

3.1 

0.3 

1.9 

6.5 

6.0 

0.9 

3.6 

17.99° (2000km) 

2.0 

2.5 

0.3 

2.2 

5.7 

5.3 

0.6 

2.8 

20° 

1.8 

2.4 

0.3 

2.4 

5.4 

5.0 

0.5 

2.7 

25° 

1.4 

1.9 

0.4 

2.7 

4.8 

4.4 

0.4 

2.7 

30° 

1.2 

1.5 

0.5 

2.9 

4.5 

4.1 

0.6 

3.1 

4o° 

1.1 

1.3 

0.6 

3.0 

4.3 

3.8 

0.9 

3.8 

50° 

1.2 

1.4 

0.5 

2.4 

4.4 

4.0 

1.0 

3.9 

60° 

1.2 

1.5 

0.3 

1.6 

4.4 

4.0 

0.9 

3.5 

70° 

1.4 

1.9 

0.3 

1.5 

4.2 

3.9 

0.8 

2.9 

80° 

1.6 

2.1 

0.3 

1.8 

3.6 

3.6 

0.5 

2.3 

90° 

1.6 

1.9 

0.3 

1.7 

3.0 

3.0 

0.4 

1.8 

100° 

1.3 

1.7 

0.3 

1.6 

2.4 

2.4 

0.3 

1.7 

110° 

1.1 

1.3 

0.2 

1.2 

2.1 

2.0 

0.4 

1.8 

120° 

1.0 

1.2 

0.2 

0.9 

2.0 

2.0 

0.5 

1.9 

130° 

1.2 

1.4 

0.2 

1.0 

2.1 

2.1 

0.5 

2.0 

i4o° 

1.3 

1.7 

0.2 

1.1 

2.0 

2.2 

0.5 

1.8 

150° 

1.2 

1.7 

0.2 

1.0 

1.7 

2.0 

0.3 

1.4 

160° 

0.9 

1.4 

0.1 

0.7 

1.1 

1.5 

0.2 

0.9 

170° 

0.3 

0.6 

0.0 

0.3 

0.4 

0.6 

0.1 

0.3 

0 

0 

00 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 
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1V9N  Nl  ID3Jj3  SWd 


Ph 

CS 

h 

a 


6o 


RADIUS  OF  NEGLECTED  ZONE 


Finally  we  consider  the  spherical  harmonics  expansion  of  Zhongolovich 


quoted  and  used  in  [Molodenskil  et  al.,  19^2^  pp„  !l45-6],  Here^ 

^2=  6l,  c^=  96,  Cj^=  12,  c^=  8,  Cg=  14,  c^=  5,  Cg=  8  mgal 

(computed  by  c  =  Ag  ^).  With  these  coefficients  we  find  values  con- 
n  n 

siderably  higher  than  those  given  above  (Table  T-^)-  They  are,  however, 
probably  too  high  since  the  development  of  Zhongolovich  is  already  out 
of  date  (it  was  computed  in  1952), 
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Table  J.k 

BMS  influence  of  the  zones  beyond  a  radius  Y  on  6  ,  d.  6^,  based  on 

O  K  M  L 

the  spherical  harmonlce  development  of  Zhongolovich,  1952  (summation 
up  to  8th  order). 


■5® 

8.99®(lOOOhm) 

10° 

13.49® (l500hm) 

''17.99®  (2000hm) 
20® 

25® 

30* 

40° 

50° 

60° 

70° 

80® 

90® 

100® 

110° 

120° 

130® 

l4o° 

150® 

160° 

170® 

180® 


Direct 


ZT 

N 

mgal 

16.2 

16.0 

14.2 

13.6 

13.8 

13.0 

12.3 

11.5 

10.7 

9.9 

10.2 

9.3 

9.1 

8.2 

8.4 

7.5 

8.1 

7.2 

8.2 

7.3 

8.2 

7.2 

7.6 

6.9 

6.8 

6.5 

5.8 

5.7 

4.8 

4.8 

4.1 

4.0 

4.0 

5.9 

4.1 

4.1 

4.0 

4.2 

3.3 

3.6 

2.0 

2.4 

0.6 

0.8 

0.0 

0.0 

Coating 


M  L 

mgal 

3.5 

13.9 

2.6 

1  r\  T 

2.4 

9.2 

1.8 

6 .6 

1.0 

4.2 

0.8 

3.7 

0.7 

4.0 

1.1 

5.3 

1.7 

7.4 

1.9 

7.5 

1.8 

6.4 

1.4 

5.0 

1.0 

4.3 

0.8 

3.7 

0.7 

3.3 

0.8 

3.4 

1.0 

3.4 

1.0 

3.7 

0.9 

3.5 

0.7 

2.7 

0.4 

1.5 

0.1 

0.5 

0.0 

0.0 
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Conclusion.  On  the  whole,  the  values  of  Fig.  la,!,  based  on  [Kaula, 
1959]  and  confirmed  by  [Uotila,  1962],  might  be  the  most  realistic 
estimates,  apart  from  the  dips.  We  should  like  to  propose  the 
smoothed  values  of  Table  7 •5. 

Table  7*5 

RMS  influence  of  the  zones  beyond  a  radius  Y  on  6  ,  5  ,  6  :  proposed 

o  W  M  L 

values . 


Direct  Coating 


Y 

=  ST 

^5"=  Ail 

0 

N 

M  L 

mgal 

N 

M  L 

mgal 

8.99°  (1000km) 

8 

8 

2 

6 

13.^9°  (1500km) 

7 

7 

2 

5 

17.99°  (2000km) 

6 

6 

2 

5 

20° 

6 

6 

2 

5 

25° 

5 

5 

2 

5 

0 

0 

ro 

5 

5 

1 

5 

45° 

5 

5 

1 

5 

60° 

h 

5 

1 

4 

90° 

3 

4 

1 

3 

120° 

2 

3 

1 

2 

150° 

2 

3 

0 

2 
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These  values  are  referring  to  a  flattening  in  the  neighborhood,  of 


1.123'^.  If  we  use  a  better  flattening,  then  these  values  should  be 
somewhat  less. 

We  see  that  for  Y  >  20°  or  30®  the  influence  of  the  distant 
o 

zones  decreases  very  slowly.  It  is,  therefore,  impractical  to  extend 

the  integration  much  farther  than  20°  (coating  method)  or  30°  (direct 

method),  unless  it  is  extended  over  the  whole  earth. 

The  effect  of  the  distant  zones  on  &  and  6.^  is  somewhat  less 

M  L 

in  the  coating  than  in  the  direct  method.  The  effect  on  6^  is  much 
less.  Eiy  using  another  method  which  will  be  described  in  the  next 
section,  however,  the  influence  of  the  distant  zones  on  6^^  can  be 
made  still  smaller. 

8.  PRACTICAL  COMPUTATION  OF  THE  GRAVITY  DISTURBANCES 
We  shall  describe  now  a  method  for  the  practical  computation  of 
gravity  disturbances.  Since  the  geoid  undulations  N  have  been  com¬ 
puted  for  large  parts  of  the  world — the  accuracy  is  as  good  as  the 
existing  gravity  material — we  can  assume  them  to  be  known.  Therefore, 
the  coating  method  can  be  applied  which  has  certain  advantages  over 
the  direct  method:  simplicity,  less  influence  of  the  distant  zones, 
etc . 

Upward  Continuation  Method.  For  the  vertical  component  6  there  is, 
-  n 

however,  a  method  which  is  still  better  than  the  coating  method. 
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Compute  the  vertical  component  at  ground  (or  at  the  geoid)  hy 

(8.1)  5  °  =Ag  +  N  . 

n  R 

St 

Since  r6  =  -r  —  is  a  harmonic  function,  Poisson's  integral  for 
n  dr 

harmonic  functions  can  be  applied  to  give 


a 

r  -  R 


r6  =  , 

n  4nR 


(*  r  R5n°  , 

JJ 


D' 

a  s 


where 


r  =  R  +  H 


D  =  (/r®  +  r^  -  2Rr  cos  Y 


Setting  as  before 


_  R  _  R 
^  "  r  ■  R  +  H  ^ 


D  =  l/l  -  2t  cos  Y  + 


we  have 


D  =  rD  . 
s 


Hence 

/f,  V  t®(l-t®)  6n° 

(8.2)  \  =  XJ  iS-8a 

o 


where  6  °  is  given  by  (8.1). 


n 


This  formula  (8.2)  is  called  "upward  continuation  integral"; 
it  can  also  be  used  for  the  upward  continuation  of  the  gravity 


anomali es  As  }  since  rAg  is  also  a  harmonic  function.  Here, 
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however,  it  is  used  for  the  upward  continuation  of  the  gravity 
disturbances. 

If  we  cOTipare  (8.2)  to  the  corresponding  formula  for  the  coating 
method, 

a 

where 

M-  =Ag  +  N  , 

then  we  find  that  the  integrand  in  (8.3)  decreases,  with  increasing 
distance  D  ,  like  1/D,  whereas  in  the  upward  continuation  integral 
(8.2)  it  decreases  much  faster,  like  I/l^  It  can,  therefore,  be 

expected  that  the  influence  of  the  distant  zones  in  the  upward  con¬ 
tinuation  method  is  much  smaller  than  in  the  coating  method. 

Indeed,  if  we  put  H  =  0,  then  6  =&°,i.e.,  the  gravity 

n  n 

disturbance  is  equal  to  the  point  value  in  P  itself  and  it  does  not 
at  all  depend  on  the  other  values  and  on  the  distant  zones.  Thus, 
for  H  =  0  the  influence  of  the  distant  zones  is  rigorously  zero  in 
the  upward  continuation  method  whereas,  as  we  have  seen  in  section  7^ 
this  is  not  at  all  the  case  in  the  coating  method  (neither  in  the  direct 
method,  as  a  matter  of  fact).  From  this  we  can  conclude  that,  even 
for  higher  elevations,  only  the  nearest  surroundings  of  P  will  be  of 
some  influence. 
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Therefore^  we  can  replace  the  terrestrial  sphere  hy  its  tangential 


plane  at  P.  If  in  this  plane  we  assume  a  rectangular  coordinate  system 
x,y  (origin  P;  the  x-axis  pointing  northward^  the  y-axis  pointing 
eastward),  we  have  instead  of  (8.2)  the  plane  formula 


(s--*)  S  J 


n  2tt 


_oo  _co 


(H®+x®+y®) 


5 — sTsTe  ^  ‘3y 


In  order  to  estimate  the  influence  of  the  distant  zones  beyond  a 

certain  distance  s  from  P  we  introduce  polar  coordinates  s,  Oi  by 
o 

X  =  s  cos  CC, 


getting 


y  =  s  sin  a, 


2tt  »  o 

.  H  .  6n 

6  =  —  7—5 — yv  3 /p  s  ds  da 

n  2tt  J  J  (H®+s®) 

0=0  s=0 


The  influence  AS  of  the  zones  beyond  the  distance  s  is  given  by 
n  o 


2n  00 


'n  2tt  ''  ''  (H®+s)^ 

0=0  s=s 

o 


77375-  s  ds  do  . 


Since  for  large  s,  [/h®  +  s®  =  s  ,  we  can  simplify  this  to 


°°  — y 

5. 


A&  =  H  r  ^  ds 

n  s* 


o 


where  5  °  is  the  average  of  8  °  over  the  circle  with  radius  s. 
n,  s  n 

If  M  is  the  upper  limit  of  the  absolute  amount  of  8  , 

n,  s 


n,  s 


<  M  . 
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we  get 


00 


/^5  <HM  f  %  =  —  M. 
n  s*  s 

s  o 

o 


To  get  some  kind  of  average  effect  AS  rather  than  this  maximum 

n 


estimate,  we  may  replace  M  by  the  average  disturbance  M  outside  the 


radius  s  which  is  certainly  close  to  Ot 


(8.5)  ET  =  —  K  . 

n  s 

o 


From  this  we  see  immediately  that,  to  get  the  same  error  for  different 


elevations,  s  must  be  proportional  to  H.  If,  e.g.,  s  =  10  H,  then 
o  o 

As"  =  0.1  M  , 
n 


Since  M  will  hardly  exceed  10  mgal,  Ab  ,  in  this  case,  will  be  smaller 

n 


than  1  mgal.  So,  it  should,  in  general,  be  sufficient  to  go  as  far  as 


10  times  the  elevation.  For  larger  H,  conditions  are  even  more 


favorable  since,  the  larger  s  is,  the  smaller  can  M  be  expected  to  be. 


Plane  formulas  analogous  to  (8.4)  also  hold  for  the  horizontal 


components  5  and  5  : 

m  JL 


(8.6) 


m 


=  -  r  r  1 
^ .i.. 


6 

m 

^Jl,° 


dx  dy  . 


But,  since  5  °  =  y?  ,  =  yT]  ,  the  use  of  these  formulas  would 

require  the  knowledge  of  the  components  T]  of  the  deflection  of  the 


vertical  in  the  entire  neighborhood  of  P  „  If  we  know  §  and  T]  ,  then 
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these  formulas  can  indeed  he  applied  with  advantage.  (To  compute  § 

and  T]  hy  differentiating  the  N-field  will  not,  in  general,  he 

accurate  enoi;igh. )  In  the  following,  however,  we  shall  regard  only 

Ag  and  n  as  known,  so  that  5  and  5  must  he  computed  hy  the 

m  a 

coating  method. 

Formulas  and  Approximations.  For  easy  reference  we  put  the  relevant 
formulas  together: 


(8. 7. a)  6^ 


4n 


CT 


do  , 


(8.7.b) 


i! 

2n 


a 


sin  Y 


^cos  a 
1  sin  a 


da 


where 


t  =  ,  D  =  /^l  -  2t  cos  Y  +  t^^ 

R+H  ' 

M.  =  Ag  +  N  . 


da  =  sin  Y  dY  dO!j 


At  least  for  the  nearest  neighborhood  of  P,  the  spherical  formulas 
can  he  approximated  hy  plane  ones; 

(8.8.a)  \  =  ^  JJ  ^  dx  dy  , 

o 


(8.8.h) 


dx  dy 
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where 


(8.9)  =  /h®  +  +  y®  . 

They  follow  from  (8.7.a,h)  by  a  series  development  with  respect 

to  Y  =  —  and  to  —  if  we  put  s  cos  a  =  x  ,  s  sin  a  =  y  . 

R  R 

In  order  to  get  an  estimate  on  the  validity  of  these  approximations; 
we  develop  1/d^  or^  which  is  equivalent, 


^  1 _ 1 _ 

"  (R+H)^D^  "(R®+r®-2Rr  cos  Y)^'® 
s 


in  this  manner,  finding 


(8.10) 

s  o 


1  3 

"  2  R  D 


®  8  R®D  ® 


“  •  •  )  . 


Here  1/D^  is  spherical  and  l/D  ^  is  the  equivalent  plane  quantity  (8.9); 
s  o 

s  a  a 

and  s  =  X  +  y 

The  first  small  quantity  in  the  bracket  is 


3  H  s‘ 


3H 


2  R  D  ^  2R  " 
o 

independent  of  the  distances.  If  H  =  63. 7  km^  this  is  0.015  or  1.5^  ; 
if  H  =  637  km  it  is  15^  .  Since  the  higher  up  we  go  the  smaller  the 
gravity  disturbances  are,  this  can  be  neglected  anyway.  The  second  term  is 


8r®d  ®  8r® 

o 

For  s  =  1000  km  this  is  O.OO3  or  0.3^  ;  for  s  =  2000  km  it  is  0.01  or 
1^  .  Thi s  can  be  neglected,  too;  the  more  so  as  1/D  ^  itself  becomes 
smaller  with  increasing  distance. 

TO 


We  see  that  up  to  20°  distance  from  P  we  can,  in  general,  use  the 
plane  formulas  (8,8,a,h)  instead  of  the  spherical  ones  (8.7.a.h).  In 
particular  this  holds  practically  always  for  the  upwaicd  continuation 
integral  (8, 8, a),  except  for  very  high  elevations  (  >  250  km,  say). 

It  is  also  useful  in  the  case  that  in  the  coating  formulas  (8.7.h)  we 
extend  the  integration  only  up  to  20°  (as  we  have  seen,  it  is  impracti¬ 
cal  to  go  much  farther  unless  the  integration  is  extended  over  the 
whole  earth,  which  of  course  is  necessary  for  highest  accuracy). 

Now  we  have  to  show  how  the  x,y  in  (8.8.a,h)  are  to  he  computed. 
The  simplest  way  is 

X  =  R(cp  -  cp  )  , 

(8.11)  ^ 


y  =  R  cos  cp  (X 
P 


another  way. 


(8.12) 


X  =  R(cp  -  cpp)  , 


y  =  R  cos  cp  (X  -  Xp)  . 

The  x's  are  the  same  in  both;  the  y’s  differ  in  that  in  (8.1l) 
we  have  the  factor  cos  cp^  and  in  (8.12),  the  factor  cos  cp  .  In  our 
latitudes  there  is,  for  large  |  cp  -  ,  a  big  difference  between  both 

formulas.  Take  cp  p=  ^0°  ^  |  cp  -  cpj  =  20°  ,  i.e.,  cp  =  20°  or  6o°.  Then 
cos  20°  =  0.91^0  , 
cos  40°  =  0.766  , 
cos  60°  =  0.500  . 
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The  differences  "between  cos  cp  and  cos  cp  are  22^  and  35^1 

P 

It  is  easily  seen  that  (8.12)  is  preferable  for  larger  distances; 
only  if  we  use  this  formula  rather  than  (8.1l)  is  the  reasoning  fol¬ 
lowing  (8.10)  applicable.  The  disadvantage  of  (8.12)  is  that  the 
5'X  5',  1°X  l°j  etc.,  blocks  are  not  rectangles  in  the  system  x,y  but 
trapezoids.  In  (8.11)  these  blocks  are  represented  by  rectangles. 
Evaluation  of  the  Integrals.  In  practice,  the  integrals  must  be  evalu¬ 
ated  by  sums.  The  formulas  (8.7.a,b)  and  8.8.a,b)  are  of  the  type 
(8. 13. a)  5  .=  J/  H  •  f(cp,X)  dq  , 

(8.13.b)  S  =  JJ  PL  .  f(x,y)  dq. 

In  5  we  have  S  °  instead  of  p;  Y,a:,D  in  (8.7.a»h)  are,  of  course, 
n  n 

functions  of  the  geographical  coordinates  cp,  X;  dq  is  the  element  of 

area,  defined  by 

dq  =  R^do  =  R^  cos  cp  dcp  dX  , 
dq  =  dx  dy, 

respectively. 

The  preceding  integrals  can  be  approximated  by  sums: 


(8.14) 

®  ^  h  \ 

where 

(8. 15. a) 
or 

\ ■  \ 

(8.15  .b) 
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where  q  is  the  area  of  a  certain  block  B  (5  X  5  or  1°  X  1°  , 
k  k 

say),  and  «p  ,  X  or  x  ,  y  are  the  coordinates  of  the  center  of 
k  k  k  k 

the  block  B 

k 

Especially  for  low  elevations  and  for  the  nearest  neighborhood  of 

P  it  is  better  to  compute  c,  in  the  integrated  form. 

k 

We  limit  ourselves  to  rectangular  coordinates  x,  y.  Then 


(Fig. 8. l) we  can  compute  c  by 


Figure  8.1 


(8.16)  c^  =  F(l)  -  F(2)  +  F(3)  -  F(4) 
where 


(8.17) 

is  the 


F(x,y)  =  f(x,y)  dx  dy 
indefinite  double  integral  of  f(x,y) 
5°F(x,y) 


f(x,y)  = 


9x9y 


so  that 
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and,  e.g.,  F(l)  is  the  value  of  F(x,y)  in  point  1. 

As  a  matter  of  fact,  this  is  equivalent  to  computing 


c^by 


(8. 18, a)  c  =  JJ  f(x,y)  dx  dy 

rather  than  by 


(8.i8.b)  C^  =  f(x^,y^)  •  q^  . 

If  f(x,y)  varies  little  in  the  block  B  considered,  then  (8.l8.a) 
and  (8.l8.b)  are  practically  equal.  If  not,  then  (8.l8.a),  i.e., 
(8.l6),  is  to  be  preferred:  if  p.  is  constant  throughout  the  block, 
then  (8.l8.a)  yields  the  correct  result  vhereas  (8.l8.b)  does  not. 
(8.l8.b)  is  in  this  case,  therefore,  subject  to  a  systematical  error. 
Let  us  nov  return  to  formulas  (8.8.a,b).  For  (8. 8. a). 


for  (8.8.b), 


f  (x,y)  = 
m 


H 

H 

2ttD 

0 

2T7(x®+y^+H®)'^'^ 

X 

X 

2ttD  ^ 
0 

"  2TT(x^+y^+H’)®'® 

y 

y 

2ttD  ® 
0 

"  2TT(x®+y®+H^)^'® 

Integration  by  (8.I7)  yields 


(8. 19. a)  F  (x,y) 


n 


tan  , 

2tt  HD 

o 
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(8.19.b)  F  (x,y)  =  -  £n  (y  +  D  )  , 
m  2tt  o 

(8.I9.C)  F^(x,y)  =  -  £n  (x  +  D^)  . 

Then,  the  coefficients  c  for  6,6,5  can  be  computed  by  (8.16). 

k  n  m  £ 

In  order  to  be  able  to  apply  (8.16),  the  figure  123^  must  be  an 
exact  rectangle.  Therefore,  x,y  must  in  this  case  be  computed  by 
(8,11)  instead  of  (8.12).  Therefore,  we  are  limited  to  the  nearest 
surroundings  of  P  where  (8.II)  and  (8.12)  are  practically  equal, 
say  to  a  rectangle  3°  ^  in  the  center  of  which  P  is  situated 
(Region  A,  p.  25).  Outside  this  rectangle,  c  can  very  well  be 

K. 

computed  by  the  simpler  formula  (8.l8.b). 

For  the  computation  of  the  contribution  of  the  innermost  zones 
it  is  necessary  to  use  as  small  blocks  as  possible.  The  smallest 

I  t 

size  is  usually  5^5  •  Even  so,  it  might  be  necessary  to  take  the 

deviation  of  the  gravity  anomalies  from  the  mean  anomaly  of  the 
respective  block  into  account. 

For  this  purpose  we  consider  the  gravity  anomaly  to  be  a  linear 
function  throughout  the  block.  Then, 

(8.20)  Z\g  =  +  aiX  +  asy  . 

Inserting  this  in  (8. 8. a)  and  (8.8.b)  and  integrating  over  one  block 

B  only  we  get 
k 
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As  =  H(a  J  +  BiJi  +  aaJp)  , 
n  o  o 


As  =  8-  J 1  ^2*^12  j 

m  o 

~  ^  +  agJss  , 

where  As  ,  As  ,  As .  are  the  contrihutions  of  the  block  B  to  5  . 

n  m  Z  k  n 


6  ,  6  and 
m  Z 

(8.21)  J  =  F(l)  -  F(2)  +  F(3)  -  F(4) 


where 


F„(x,y)  .  -  ;/ 


dx  dy  1 

D  ' 
o 


— g—  =  -  tan 

D  ^  2ttH 


(-1)  2^  _  i  F  y) 

HD  K  ' 

o 


Fi(x,y)  =  ^  n 


X  dx  dy 
D  ^ 


—  ^n(y+D  )  =  F  (x,y)  , 
2tt  o  m 


(8.22) 


F2(x,y)  =  ^  a 


o 


Fii(x,y)  =  ^  ;;  ^  '  h  in(x+D  )  -  f-  tan 


H  ^  (-1)  xy 


2tt 


2tt 


2tt 


HD  ' 
o 


Fi2(x,y)=^  JJ* 


xy  dx  dy  Dp 
D  ®  "  "  2n 


/  \  1  f'py^dxdy  x^/^^  H,  (-1)  xy 

Fpj>(x,y)  =  —  - 3 -  =  —  <Cn(y+D  )  -  —  tan  . 

J  2tt  D  2tt  o  2tt  HD 


Consider  now  the  contribution  of  the  innermost  block.  Let  the 
sides  of  the  block  be  2a,  2b,  so  that  the  points  1,  2,  3^  ^  have  the 


coordinates : 


l(-a,  -b) 
3(+a,  +b) 


2  (-a,  +b) 
i+(+a,  -b)  . 
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Then, 

Ji  = 

Ja  =  <7i2 

and 

(8. 23. a) 

Ab 

=  a  HJ  , 

n 

o  o 

(8.23.b) 

Ab 

m 

=  > 

(8.23.c) 

=  ^■alaa  > 

where 

J  ,2.  (.1)^ 

o  nH  HD 

0,1 


(8,24) 


■*ss 


b  . 

Dn  1 

a 

2H 

=  —  in  — 

—  - 

TT 

^0  >1 

• 

a 

n 

=  ^  in  ^ 

+ 

_b 

TT 

^0,1 

- 

b 

TT 

tan 


HDo ,  1 


HDo  1 


Do,i  =  |/a*  +  D®  +  H®  , 


We  see  that  for  the  central  block  A5  is  the  same  as  if  we 


n 

would  use  formula  (8.l6).  For  ^6  and  As  (8.l6)  yields  zero 

m  n 

whereas  the  correct  values  are  given  by  (8.23.b,c). 

In  order  to  estimate  numerically  the  contribution  of  the 

innermost  compartment,  Ab  and  Ab  ,  we  take  it  as  a  square  with 

m  Z 

sides  s  ,  so  that 


a  =  b  =  s  ;  Do,i  =  / 2s*  +  H*  , 


Since,  then,  =  Jgs  it  is  sufficient  to  consider  A b^  only. 
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We  have 


(8.25) 


Jii 


Pq  1 1  +  s 
Do , 1  -  s 


2H  (-1)  s® 

—  tan  - 

n  HDo, 1 


we  take 


s  =  4-  km 

I  I 

which  corresponds  to  the  average  size  of  5  k  5  blocks, 
gives  the  values  of  for  several  elevations  H  . 

Table  8,1 

Factor  for  s  =  4  km  as  a  function  of  elevation 

H  Jii 

km  km 


0 

2.24 

1 

1.41 

2 

0.87 

5 

0.23 

10 

0.04 

20 

0.01 

50 

0.00 

100 

0.00 

If  we  have  a  north- south  gradient  ai  =  1  mgal  per  1 
Jii  is  numerically  equal  to  the  effect  A5^  in  mgal;  e.g. 


Table  8.1 

H 


km,  then 
,  for  H  =  0 


we  have 


Since  also  gradients  of  2  or  even  3  mgal/km  are  possible,  the 

effect  of  the  gradient  on  the  horizontal  components  A5  and  Ab 

m  Z 

can  be  considerable  for  lower  elevations,  H  <  5  km.  For  H  >  5  km, 
the  effect  of  the  gradient  is  negligible. 

We  have  also  to  investigate  the  influence  of  the  eight  blocks 
adjacent  to  the  central  one.  Since  the  effect  is  a  maximum  for 
H  =  0  ,  we  consider  this  case  only.  We  again  take  the  block  to  be  a 
square,  but  with  arbitrary  sides.  We  find  that  the  effect  of  the 
innermost  block  is 

(8.26)  As  =  0.561  aiS 

m 

and  the  combined  effect  of  the  eight  surrounding  blocks  is  only 

(8.27)  AI6  =  -  0.071  ais  . 

in 

I  I 

If  the  gradient  ai  =  1  mgal/km  and  s  =  4  km  (5  k  5  blocks),  then 

A'5  =-0.3  mgal  . 

m 

These  values  hold  for  H  =  0  ,  but  even  here  the  effect  of  the  gradient 
in  the  surrounding  blocks  will  be  negligible  in  most  cases,  even  if 
it  must  be  taken  into  account  for  the  central  block  (for  H  <  5  km). 

These  considerations  are  equally  valid  for  the  component  6^  ; 
we  have  only  to  replace  ai  by  the  east-west  gradient  ag  .  For  the 
vertical  component  6  ,  the  influence  of  the  gradients  ai  and  ag 

will  always  be  negligible. 
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An  easy  way  of  computing  the  gradients  aj.  and  wMch  can  also 
be  used  for  high-speed  computations  is  (Fig.  8.2): 


Fj-gure  8.2 

Station  and  Substations 


(8.28) 


/\gi  through  Ag4  are  the  gravity  anomalies  in  the  centers  of  the 
adjacent  blocks  1  to  4,  which  we  call  substations.  These  point 


anomalies  can  also  be  replaced  by  the  mean  anomalies  of  the  respective 
blocks,  provided  the  gravity  anomalies  are  sufficiently  linear  (this 
is  presupposed  in  (8.28)). 


Station  and  Substations.  For  reasons  of  symmetry  the  station  P 

at  which  6,6,6  are  computed  must  be  at  the  center  of  the 
n  m  £ 

Y  I 

central  block.  If  we,  however,  use  given  5  X  5  mean  anomalies 


for  which  the  division  into  blocks  is  fixed  a  priori,  it  is  best  to 
select  5  substations  (centers  of  adjacent  blocks)  0,  1,  2,  3,  ^  (Fig.  8.2) 
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and  compute  the  gravity  disturbance  in  those  three  substations  which 
are  the  corners  of  a  triangle  containing  P  in  its  interior  (O,  1, 

4  in  Fig.  8.2).'  The  value  for  P  is  then  interpolated  between  thesfe 
three  substations.  This  procedure  is  particularly  necessary  for 
elevations  smaller  than  about  70  km. 

Simimary  of  Formulas  and  Outline  of  Computation.  In  the  following 
we  shall  give  a  computational  method  valid  for  any  elevation  from 
0  to  several  hundred  kilometers.  It  is,  therefore,  suitable  for  the 
computation  of  the  gravity ■ vector  along  rocket  trajectories. 

According  top. 25  we  use  the  following  regions 

A  9i  =  .3°  ,  '  , 

5  X  5  blocks, 

X2  -  Xi  =  4° 

B  %  -  9i  =  7°  ,  , 

outside  A,  20  X  20  blocks, 

Xa  -  Xi  =  9° 

C  CPg  -  CPi  =  25 

outside  B,  1°  X  1°  blocks, 

Xs  -  Xi  =  30°  .  .  .... 

D  outside  C,  5°  X  5°  blocks. 


1 ,  For  each  block  compute 


|i  =  Ag  +  0.231  N  ,  5“  =  Ag  +  0.308  N  ; 

Ag:  mean  free-air  anomalies,  N:  geoid  undulations;  As^ 
mgal,  N  in  meters. 

2,  For  each  station  P  select  3  substations,  according  to  Fig.  8.2. 
One  substation,  0,  is  the  center  of  the  block  in  which  P  lies,  the 
other  two  are  centers  of  adjacent  blocks  so  that  the  three  substations 
form  the  comers  of  a  triangle  containing  P  .  The  elevation  of  each 
substation  is  equal  to  the  elevation  H  of  P  . 

The  gravity  disturbances  are  computed  at  each  substation  and 
are  interpolated  for  P  . 

In  the  following,  only  the  substations  will  be  considered;  they 
will  be  denoted  by  Q  . 

3.  Now  we  consider  one  definite  substation  Q  only. 

Zone  A. 

Compute  for  all  grid  points  (  =  corners  of  blocks): 

X  =  R(cp  -  «p^)  , 
y  =  R  cos  cp^  (X  -  X^)  ; 

D  =  /  X*  +  y*  +  ; 

o  ' 
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F 

n 


(-1)  ^ 

hD 


} 

o 


F 

m 


(y  +  D) 

O 


F  =  -  —  in  (x  +  D)  . 
i  2n  '  o 


Compute  for  each  block 


c  =  F  (1)  -  F  (2)  +  F  (3)  -  F  (4)  , 

n  n  n  n  n 


3)  -  F  (4)  , 

mm  m  m  m 


c,  =  F  ^1)  -  F,(2)  +  F  (3)  -  F  (4)  . 

i  V  ^  '  ‘C  ’  V  ‘ 

Be  oure  that  the  grid  points  1,  2,  3^  ^  are  arranged  as  in  Fig.  8,1. 

For  elevations  H  >  40  km  the  coefficients  c  can  also  be  com¬ 
puted  as  in  zone  B. 

Zone  B 

Compute  for  all  block  centers 

X  =  R(cp  -  cp^)  ^ 

y  =  R  cos  cp(X  -  \  )  ; 

D  =  /  x^  +  y^  +  H®  ; 
o 


H 


c 


n  2nDl 


m  2nDr 


c  ,  = 


'I  2ttD^ 
o 


q  is  the  area  of  the  block. 
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Zone  C  can  be  treated  either  like  Zone  B  or  like  Zone  D. 


Zone  D 

If  many  points  are  to  be  computed  then  the  following  computations  need 
not  be  performed  for  all  points.  In  this  case  it  is  sufficient  to 
compute  the  effect  of  2ione  D  for  a  few  points  only  and  interpolate 
between  these, 

CcHupute  ; 


Y  =  cos  ^^[sln  cp^  sin  (p  +  cos  <p^  cos  cp  cos  (X  -  X^)]  , 

Q  <P  -  cos  tP  cos  (X  ■»  Xq)  j 

cos  (p  sin  (X  -  Xq)  " 


D  =  ^  1  -  2t  cos  Y  +  t*  ; 


c 

n 


H  ^  200  km 


q  ,  H  >  200  km 


sin  ¥  cos  a  » 


<1  » 


t® 


4.  For  the  substation  Q  considered  compute  by  suawation  over  all 
blocks  1  *. 
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S  =  E  c  5°  , 

n  ^  n,i  n,i 

5  =  E  c  .  M-. , 

m  .  m,i  1 


5  ,  =  ^  c  .  u.  . 

For  elevations  <  5  km  we  have  to  add,  to  5  and  to  &  ,  the  terms 

m  Z 

^  and  respectively: 

m  Z 


As  =  ai 
m 


r  h  ^  Di+a 
[  -  Zn  — — 
rr  Di-a 


2H  ^  (-1)  ah 

—  tan  —  J , 

TT  ■  HDi  ' 


AS^  =  ag 


r  a  .  Di+h 
[  -  Zn 

TJ  Dj-D 


2H  (-1)  ah  ■ 
—  tan  —  J , 

TT  HDi  ^ 


where  2a, 2h  are  the  sides  of  the  central  5'  k  5 '  block  and 


Di  =  /a^+  h^+ 

Agi-  Aga  Ags-  Ag4 

=  4a  ^  =  4h 


These  formulas  for  a^  and  ag  are  valid  for  the  central  substation  0 
(Fig.  8.2);  for  the  other  substations  they  are  to  be  modified  in  an 
evident  way. 


5.  Now  denoting  by  x,y,z  the  coordinates  in  the  geocentric  system 
defined  on  p.  2  (z-axis  =  axis  of  rotation  of  the  earth)  compute  the 
components  of  the  gravity  disturbance  in  this  system  by 
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8  =  -  cos 

X 

«P 

cos 

X  •  8  -  sin 

n 

9 

cos 

\  • 

8  -  sin 

m 

CO 

• 

8  =  -  cos 

y 

9 

sin 

X  •  8  -  sin 

n 

9 

sin 

X  • 

8  +  cos 

m 

CO 

• 

5  =  -  sin  cp  •  8  +  cos  cp  •  8  . 

z  n  m 

These  components  8,8,8  are  positive  away  from  the  earth' s  center; 
X  y  z  - 

cp,X  refer  to  the  substation  Q  . 

6.  This  procedure  must  he  repeated  for  all  substations  and  the  values 

of  8  ,  8  ,  8  for  the  station  P  are  interpolated. 

X  y  z 

For  elevations  H  >  70  km,  it  is  in  general  sufficient  to  use  the 
central  substation  0  only  and  to  put  the  values  at  P  equal  to  the 
values  at  0  . 


7.  Computation  of  the  components  of  the  gravity  vector; 


g  =  ty  -  8  , 
X  X  X 


g  =  y  -  5  , 
y  y  y 


g  =  y  -  5  , 

z  z  z 


y  y  y  i  y  positive  toward  the  earth’s  center. 

X  y  z  X  y  z  ■  -  ■  ■■ 
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